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Abstract 


Geolocation  is  a  common  application  for  satellite  systems.  This  involves  estimating 
an  object’s  location  (herein  called  the  subject)  based  on  noisy  satellite  data.  Many 
geolocation  methods  exist;  however,  none  are  tailored  specifically  for  the  unique  problems 
faced  by  satellite  systems.  Some  satellites  are  so  far  from  the  subject  being  localized 
that  by  the  time  the  satellite  receives  a  signal  from  the  subject  it  might  have  moved 
appreciably.  Furthermore,  some  satellites  or  terrestial  sensors  may  be  much  closer  to 
the  subject  than  others.  Therefore,  sensors  may  need  to  be  weighted  based  upon  their 
distance  to  the  subject  being  localized.  In  addition,  even  if  a  subject  can  be  localized, 
the  confidence  in  this  localization  may  be  unknown.  Non-linear  optimization  is  proposed, 
implemented,  and  analyzed  as  a  means  of  geolocating  objects  and  providing  confidence 
estimates  from  passive  satellite  line-of-sight  data.  Non-linear  optimization  requires  an 
initial  estimate.  This  estimate  is  provided  by  a  triangulation  method.  The  non-linear 
optimization  then  improves  upon  this  estimate  iteratively  by  finding  estimates  that  are  more 
likely  to  have  produced  the  observed  line-of-sight  measurements.  The  covariance  matrix  of 
the  geolocation  parameters  being  estimated  is  naturally  produced  by  the  optimization  which 
provides  quantified  confidence  in  the  geolocation  estimate.  Simulations  are  developed  to 
provide  a  means  of  evaluating  the  performance  of  the  non-linear  optimization  algorithm. 
It  was  found  that  non-linear  optimization  can  reduce  the  average  error  in  geolocation 
estimates,  provide  improved  estimation  confidence,  and  accurately  estimate  its  geolocation 
confidence  for  some  subjects.  The  results  from  the  theoretical  development  of  the  non¬ 
linear  optimization  algorithm  and  its  simulated  performance  is  quantified  and  discussed. 
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NON-LINEAR  OPTIMIZATION  APPLIED  TO  ANGLE- OF- ARRIVAL 


SATELLITE-BASED  GEOLOCATION 

I.  Introduction 

THis  chapter  introduces  the  motivation  behind  this  research.  It  describes  the  problem 
at  hand  and  the  associated  research  question.  It  discusses  the  assumptions  and 
limitations  to  this  research.  Lastly,  it  provides  an  outline  for  the  rest  of  the  thesis. 

In  this  research,  Non-Linear  Optimization  (NLO)  algorithms  are  developed  and 
investigated  as  a  means  of  localizing  an  object  from  Angle-of-Arrival  (AOA)  data. 
Localizing  an  object  means  to  estimate  its  position  in  space.  For  the  purposes  of  this 
research,  the  object  which  is  being  localized  will  be  referred  to  as  the  subject.  While  AOA 
data  from  any  source  may  be  used,  this  research  focuses  on  data  passively  obtained  by 
satellites.  Thus,  the  sensors  do  not  transmit  and  receive  an  echo  from  the  subject.  Rather, 
the  sensor  observes  an  emission  from  the  subject.  Localization  provides  an  estimate  of  the 
subject’s  location,  but  without  an  estimate  of  the  confidence  in  the  subject  location,  this 
information  may  have  limited  utility. 

Therefore,  NLO  is  also  investigated  for  its  ability  to  estimate  the  confidence  it 
provides  in  its  localization.  Furthermore,  the  algorithm  must  provide  a  means  of  intuitively 
conveying  this  information  to  a  user  via  some  visualization. 

To  conduct  this  research,  the  NLO’s  ability  to  localize  an  object  and  represent  its  error 
is  theoretically  developed.  The  theoretical  performance  is  compared  to  its  performance 
in  simulated  scenarios.  A  Scenario  consists  of  a  subject,  observing  sensors,  and  their 
measurements.  Two  types  of  scenarios  are  produced:  scenarios  for  verification  and 
performance  evaluation. 


1 


The  first  type  of  simulation  is  used  to  verify  that  the  NLO’s  theoretical  performance 
matches  its  simulated  performance.  NLO  is  based  on  several  assumptions.  Therefore, 
in  scenarios  where  these  assumptions  are  true,  the  theoretical  results  and  simulated 
results  should  match.  This  match  verifies  that  the  NLO  has  been  correctly  implemented. 
Furthermore,  these  types  of  simulations  are  useful  for  comparing  the  NLOs.  The  NLOs 
include  different  assumptions  which  add  to  their  complexity.  Thus,  these  results  may 
illustrate  how  much  each  assumption  separately  impacts  their  ability  to  localize  the  subject 
and  represent  the  confidence  in  their  geolocation  estimates. 

The  second  type  of  simulation  models  a  scenario  where  the  NLOs’  assumptions 
are  imperfect.  Therefore,  these  simulations  are  used  to  evaluate  the  performances  of 
NLOs  in  scenarios  where  their  assumptions  are  approximations.  The  NLOs  will  produce 
less  accurate  estimates  if  their  assumptions  are  not  perfectly  true.  By  comparing  the 
performances  of  the  NLOs  in  this  scenario,  the  error  resulting  from  unconsidered  factors 
may  be  investigated. 

This  research  provides  the  theoretical  development  for  three  NLO  algorithms  of 
varying  complexity  for  localization  using  AOA  data.  It  describes  how  these  algorithms 
are  verified  and  how  their  performances  are  evaluated.  The  performance  and  development 
of  the  confidence  visualization  scheme  are  also  given.  The  next  section  describes  how  this 
research  started  and  why  it  is  a  relevant  effort. 

1.1  Background 

In  many  applications,  it  is  important  to  estimate  the  positions  of  objects  and  provide 
the  accuracy  in  these  estimates  via  a  process  commonly  referred  to  as  geolocation  or 
localization  [23],  [27],  [34].  Geolocation  provides  information  on  objects  world-wide.  One 
of  the  most  common  means  of  estimating  the  position  of  an  object  is  via  the  use  of  a  Global 
Navigation  Satellite  System  (GNSS)  such  as  the  Global  Positioning  System  (GPS)  [16]. 
Each  GPS  satellite  transmits  a  signal  that  includes  the  time  the  signal  was  sent,  along  with 
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the  satellite’s  position  at  that  time.  A  receiver  may  then  use  this  information  from  multiple 
satellites  to  produce  an  estimate  of  the  object’s  position  by  a  method  known  as  trilateration 
[26], 

There  is  an  inverse  problem  of  sorts  known  as  passive  geolocation  or  source 
localization  [10].  In  this  problem,  the  emitter’s  position  is  unknown  and  must  be  estimated 
using  measurements  from  multiple  receivers.  It  is  often  assumed  that  the  emitter  does  not 
transmit  its  location.  Such  an  emitter  is  said  to  be  non-cooperative.  In  this  research,  the 
subject  is  a  non-cooperative  emitter.  Because  the  subject  is  non-cooperative,  the  receivers 
must  produce  the  information  needed  to  localize  the  source. 

Therefore,  differences  in  the  signals  received  by  multiple  observing  sensors  must  be 
exploited  for  source  localization.  The  two  kinds  of  data  are  used  predominantly  in  passive 
geolocation:  Time  Difference  of  Arrival  (TDOA)  and  AOA  (or  equivalently  Direction  of 
Arrival  (DOA))  measurements  [4],  [11].  TDOA  techniques  use  the  differences  between  the 
times  when  an  emitter’s  signal  is  received  by  multiple  sensors  to  perform  localization  [29]. 
AOA  methods  use  the  direction  or  Line-of-Sight  (LOS)  to  the  emitter.  If  these  LOS  are 
extended,  they  should  meet  at  a  point.  Therefore,  the  point  nearest  to  the  intersection  of  the 
LOS  may  be  used  to  estimate  the  emitter’s  location. 

Passive  geolocation  is  useful  in  many  applications.  For  example,  it  may  be  useful 
where  GPS  is  not  available.  It  also  provides  a  method  for  geolocating  non-cooperative 
subjects  [29],  [32].  The  ability  to  localize  non-cooperative  subjects  has  led  to  research 
in  satellite-based  passive  geolocation  systems  that  use  LOS  information  [17],  [33].  LOS 
geolocation  must  be  used  in  some  cases  where  the  TDOA  cannot  be  determined.  For 
example,  if  a  signal  is  emitted  continuously,  then  the  start  of  the  signal  may  be  unknown. 
The  phase  of  the  signal  could  also  be  exploited  to  determine  the  TDOA;  however,  if  the 
frequency  of  the  signal  is  too  high,  then  TDOA  methods  may  be  impossible.  Therefore, 
AOA  geolocation  techniques  may  provide  geolocation  when  other  methods  cannot  be  used. 
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There  are  several  standard  LOS  geolocation  schemes.  These  often  employ  some  form 
of  triangulation  [20],  [29].  While  these  schemes  are  theoretically  and  computationally 
simple,  they  do  not  account  for  the  unique  challenges  faced  by  satellite  geolocation 
systems.  These  challenges  are  discussed  below. 

Challenge  1:  the  further  a  sensor  is  from  the  subject,  the  more  an  error  in  its  LOS 
produces  an  error  in  geolocation.  This  fact  is  illustrated  in  Figure  1.1.  Errors  in  the  LOS 
measured  by  sensors  that  are  close  to  the  subject  produce  smaller  geolocation  errors  than 
sensors  that  are  further  from  the  subject.  This  issue  is  particularly  relevant  for  satellite- 
based  sensors  where  the  distances  involved  may  be  large.  Thus,  the  closer  sensor  should  be 
weighted  more  heavily  than  a  far-distant  satellite  since  its  error  propagates  over  a  shorter 
distance.  Because  the  weighting  for  each  sensor  depends  on  its  accuracy  and  the  emitter’s 
location  (which  is  unknown),  the  weights  cannot  be  calculated. 


Figure  1.1:  Error  propagation  and  triangulation.  Both  satellites  have  the  same  angular 
error.  However,  the  satellite  that  is  further  from  the  subject  (shown  in  red)  has  this  error 
propagated  over  a  larger  distance. 


Challenge  2:  Another  factor  ignored  by  triangulation  is  that  satellites  may  be  so  far 
from  the  subject  that  the  subject  may  have  moved  a  significant  distance  by  the  time  its 
emission  is  detected  by  the  satellites.  This  problem  is  referred  to  as  the  time  delay,  and  it 
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is  illustrated  in  Figure  1.2.  Therefore,  the  satellite’s  LOS  is  a  measure  of  where  the  emitter 
used  to  be  rather  than  where  it  is  when  the  sensor  produced  its  measurement. 


Figure  1.2:  The  time  delay  problem.  The  left-hand  frame  shows  the  subjecft  when  it 
emitted  a  signal  (shown  in  yellow).  By  the  time  this  signal  is  received  by  the  sensor  (shown 
in  frame  2),  the  subject  has  moved  a  significant  distance. 


Challenge  3:  Furthermore,  satellites  may  not  produce  LOS  measurements  simultane¬ 
ously.  This  problem  is  illustrated  in  Figure  1.3.  Finding  the  point  nearest  the  intersection  of 
these  LOS  assumes  that  the  subject  is  in  one  place  when  these  measurement  are  produced. 
Each  satellite’s  LOS  is  measured  at  different  times.  For  a  moving  subject,  these  LOS  mea¬ 
surements  point  toward  the  subject  at  different  points  along  its  path.  Therefore,  the  subject 
must  be  localized  at  a  chosen  time.  Thus  measurements  produced  closer  to  this  time  should 
be  given  more  weight. 

Challenge  4:  Furthermore,  triangulation  may  not  calculate  the  confidence  in  its 
geolocation.  Therefore,  there  is  no  perfect  measure  of  how  accurate  or  precise  the 
localization  may  be.  For  these  reasons,  it  is  desirable  to  find  another  geolocation  method 
which  resolves  the  issues  associated  with  satellites. 

It  has  been  noted  [24]  that  geolocation  typically  has  a  nonlinear  nature.  Thus,  the 
NLO  has  been  considered  for  geolocation  [9] .  However,  NLO  has  never  been  applied  to 
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Figure  1.3:  Triangulation  problem  of  asynchronous  LOS  measurement  times.  The  LOS  to 
the  subject  has  been  measured  at  three  separate  points  along  its  track  (the  red  line). 


geolocating  satellite  systems  using  LOS  data.  In  addition,  NLO’s  ability  to  resolve  the 
aforementioned  shortcomings  has  not  been  investigated. 

1.1.1  Research  Question  &  Significance. 

The  focus  of  this  research  is  to  investigate  two  queries:  First,  how  well  does  NLO 
geolocate  the  subject  in  a  satellite  system  that  measures  AOA?  Several  NLOs  might  be  used 
(three  are  considered  here).  This  question  then  extends  to  asking  which  NLO(s)  perform 
the  best.  Secondly,  is  there  a  useful  means  of  estimating  and  visualizing  the  confidence  in 
the  geolocation  estimates  it  provides?  If  the  NLO  is  capable  of  overcoming  triangulation’s 
shortcomings,  it  may  be  able  to  provide  users  with  several  benefits. 

First,  if  NLO  compensates  for  the  sources  of  error  in  satellite  geolocation,  then 
position  estimates  will  be  more  accurate.  The  NLO  may  also  provide  better  confidence 
than  triangulation  in  its  estimate,  thereby  providing  a  benefit  to  the  user.  This  improved 
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confidence  may  allow  for  other,  more  accurate,  sensors  to  hone  in  on  the  emitter  more 
quickly. 

1.2  Research  Assumptions  and  Domain 

This  research  involves  multiple  assumptions  and  practical  limitations.  In  addition,  this 
research  evaluates  the  research  question  only  for  a  specific  set  of  scenarios.  The  following 
section  discusses  the  research  assumptions. 

1.2.1  Assumptions. 

For  most  passive  geolocation  techniques,  knowledge  of  the  sensors’  positions 
(satellites  for  the  purposes  of  this  research)  is  requisite.  In  reality,  there  may  be  some 
uncertainty  in  sensor  locations.  For  the  purposes  of  this  research,  this  issue  is  ignored  and 
it  is  assumed  that  the  satellites’  positions  are  known  absolutely. 

The  Lines-of-Sight  (LOS)  to  the  subject  is  also  uncertain.  Measurements  typically 
have  two  sorts  of  errors:  systematic  and  stochastic  errors.  In  this  research,  it  is  assumed 
that  the  errors  are  stochastic  and  uncorrelated.  Stochastic  errors  are  typically  described 
by  a  probability  distribution  function.  A  Gaussian  distribution  is  used  as  a  common  first 
assumption.  This  distribution  is  assumed  for  random  variables  in  this  research. 

The  distribution  is  described  in  mathematical  detail  in  Chapter  2.  It  is  also  assumed 
that  these  measurements  are  Independent  and  Identically  Distributed  (IID).  That  is,  each 
measurement  is  independent  of  the  measurements  which  came  before  and  after  it.  Gaussian 
distributed  random  variables  are  often  assumed  to  be  IID. 

In  addition  to  assumptions  about  the  measurements  and  the  sensors,  assumptions  are 
made  about  the  subject.  Geolocation  performance  may  vary  with  the  type  of  subject.  To 
simplify  this  research,  only  one  type  of  subject  is  considered.  Given  the  discussion  on  the 
sources  of  error  discussed  in  Section  1.1,  it  is  assumed  that  a  Space  Vehicle  (SV)  is  the 
most  challenging  subject  on  average.  An  SV  is  most  challenging  in  the  sense  that  SVs  are 
assumed  to  result  in  worse  accuracy  and  confidence  than  any  other  subject  on  average.  SVs 
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produce  the  worst  geolocation  performance  because  they  are  both  the  fastest  subjects  that 
might  be  observed  and  their  flight  paths  are  nonlinear. 

It  is  also  assumed  that  the  times  at  which  each  satellite  receives  a  signal  is  exactly 
known.  This  assumption  is  justified  for  the  following  reason:  the  times  at  which  the 
satellites  produce  their  measurements  matter  only  because  the  emitter  may  be  in  motion. 
Therefore,  if  the  reception  time  is  off  by  some  amount,  then  the  emitter  may  have  moved 
some  distance  during  that  time.  Given  the  accuracy  of  satellite  clocks,  even  the  fastest 
emitter  could  not  move  an  appreciable  distance  over  the  duration  of  the  timing  error.  In 
addition  to  these  assumptions,  there  are  practical  constraints  and  limitations. 

1.2.2  Limitations. 

NLO  is  fundamentally  limited  in  that  it  optimizes  over  a  model  of  the  scenario.  The 
model  describes  the  mathematical  relationship  between  the  subject’s  motion  and  satellite 
positions  and  the  LOS  measurements.  The  model  will  not  perfectly  represent  the  subject’s 
motion. 

Real  subjects  will  have  velocity,  acceleration,  jerk,  jounce,  and  other  higher  order 
moments.  For  example,  if  the  NLO  is  constructed  to  account  for  velocity,  then  it  assumes 
that  the  emitter  is  approximately  moving  in  a  straight  line  over  the  duration  of  the 
measurements  used  for  geolocation.  Because  there  is  a  limit  to  the  number  of  parameters 
that  can  be  included,  the  NLO  cannot  include  arbitrary  subject  motion  in  its  model.  In  this 
research,  only  the  position  and  velocity  of  the  emitter  and  the  time  delay  from  the  emitter 
to  the  satellite  are  optimized.  This  thesis  evaluates  whether  satisfactory  results  may  be 
achieved  given  these  NLO’s  models. 

A  further  practical  limitation  of  this  research  is  real  LOS  data.  No  real  data  was 
available  for  this  research.  Likewise,  typical  confidence  in  LOS  is  unknown.  Therefore,  the 
NLO  is  evaluated  on  simulated  data  with  an  arbitrary  confidence  in  the  LOS  measurements. 
The  scope  is  this  research  is  described  in  the  next  section. 
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1.2.3  Scope. 

This  research  is  strictly  focused  on  satellite  based  geolocation.  Other  platforms 
could  be  used  in  conjunction  with  a  satellite  system.  The  effects  of  terrestrial  sensors 
on  geolocation  performance  are  not  considered.  Furthermore,  TDOA  and  LOS  information 
could  be  available.  In  this  case  both  could  be  used  in  conjunction  to  provide  geolocation 
estimates.  This  is  not  considered. 

Neither  is  the  full  range  of  possible  NLO  algorithms.  Only  three  NLO  algorithms 
are  considered  here.  It  is  anticipated  that  these  three  optimize  the  most  significant  unique 
sources  of  error  as  previously  described. 

Another  geolocation  task  that  this  research  ignores  is  distinguishing  LOS  information 
for  multiple  subjects.  It  may  be  ambiguous  which  LOS  measurements  relate  to  which 
subject,  and  this  difficulty  is  compounded  if  the  number  of  emitters  is  unknown.  For 
example,  two  emitters  could  be  at  nearly  the  same  LOS  relative  to  a  satellite,  but  their 
distances  to  the  satellite  could  vary  drastically.  It  would  be  difficult  to  determine  which 
emitter  corresponds  to  which  LOS. 

The  scope  of  this  research  also  excludes  the  computational  speed  of  the  NLO.  There 
may  be  requirements  regarding  the  speed  at  which  a  geolocation  algorithm  performs.  The 
speed  of  the  NLO  is  recorded,  but  it  is  not  used  as  grounds  for  evaluating  the  NLO’s 
performance. 

1.3  Evaluation  Methodology 

To  properly  evaluate  the  benefits  provided  by  the  NLO,  its  geolocation  performance  is 
evaluated  in  four  ways:  the  mean  absolute  error  in  the  geolocation  estimate,  the  bias  in  its 
estimate,  the  confidence  in  its  estimate,  and  the  accuracy  of  its  confidence  estimate.  These 
are  collectively  referred  to  as  performance  of  the  geolocation  algorithm.  These  values  are 
anticipated  theoretically  and  calculated  in  simulated  scenarios. 


9 


The  theoretical  properties  of  the  NLO’s  performance  metrics  are  developed  analyt¬ 
ically.  The  theoretical  performance  is  then  verified  by  simulated  scenarios  in  which  the 
NLO’s  assumptions  are  exact.  Upon  verification  of  the  NLOs  theoretical  performance,  the 
NLO  is  used  in  simulated  scenarios  where  the  NLO’s  assumptions  are  approximately  true. 
The  NLO’s  performances  in  these  scenarios  are  compared  against  each  other  and  another 
common  LOS  geolocation  algorithm:  a  Triangulation  Algorithm  (TA).  The  TA  does  not 
take  into  account  the  unique  sources  of  error  posed  by  satellite  systems.  Thus,  the  com¬ 
parison  between  the  TA  and  NLO  provides  data  on  the  performance  increase  or  decrease 
provided  by  optimizing  for  these  sources  of  error.  The  comparison  between  the  NLO’s  per¬ 
formances  provides  the  geolocation  improvement  resulting  from  adding  model  complexity 
to  the  NLO. 

To  estimate  the  mean  absolute  error  and  bias  of  the  algorithms,  truth  data  and  noisy 
LOS  measurements  are  simulated.  The  truth  data  are  the  true  LOS  to  the  subject  and  the 
subject’s  position  (and  velocity  if  relevant).  The  TA  and  NLO  use  the  noisy  data,  and  their 
geolocation  estimates  are  compared  to  the  truth  data.  This  comparison  provides  a  metric 
for  quantifying  the  accuracy  provided  by  the  NLO  and  TA. 

The  confidences  provided  by  the  NLOs  and  TA  are  described  by  the  variance  in  their 
estimates.  The  comparison  between  their  true  confidence  and  the  confidence  they  estimate 
is  also  calculated.  This  is  the  accuracy  in  their  confidence  estimates.  These  metrics 
are  calculated  from  synthesized  scenarios.  Monte  Carlo  Simulations  (MCS)  are  used  to 
provide  the  algorithms’  true  confidences.  The  NLOs  have  a  means  of  calculating  their 
confidence  without  MCS.  Therefore,  their  confidence  estimates  are  compared  against  the 
confidence  estimates  produced  from  MCS.  The  confidence  and  accuracy  of  the  confidence 
estimates  provided  by  the  NLOs  are  compared.  As  mentioned  in  Section  1.2,  three  NLOs 
of  various  complexity  are  evaluated.  Therefore,  the  benefit  of  adding  more  parameters  (and 
complexity)  to  the  NLO  is  analyzed. 
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To  find  the  performance  differences  provided  by  the  NLOs,  the  NLO  is  developed  in 
three  stages  of  optimization.  The  first  NLO  optimizes  the  emitter’s  position,  the  second 
NLO  optimizes  position  and  velocity,  and  the  final  NLO  includes  the  time  delay.  This 
research  compares  the  accuracy  and  confidence  provided  by  these  three  NLOs  to  evaluate 
the  performance  difference  resulting  from  the  added  complexity. 

1.4  Overview 

Chapter  2  begins  with  an  explanation  of  the  data  used  by  the  geolocation  algorithms. 
It  proceeds  by  explaining  how  simulations  create  these  data  and  how  the  most  challenging 
subject,  the  SV  is  synthesized.  The  chapter  proceeds  by  describing  and  deriving  the 
mathematics  behind  the  TA  and  general  NLO  algorithms.  The  chapter  closes  by  explaining 
how  the  confidence  estimates  provided  by  the  geolocation  algorithms  are  calculated  and 
visualized. 

Chapter  3  discusses  how  the  theoretical  building  blocks  described  in  Chapter  2  come 
together  in  simulations.  It  describes  the  combined  synthesis  of  the  satellites,  their  LOS 
measurements,  and  the  subject.  Next,  it  is  explained  how  specific  types  of  scenarios  are 
simulated  to  verify  the  functionality  of  the  NLOs.  This  is  followed  by  an  explanation  of 
each  NLO  and  the  subtleties  of  their  operation.  The  chapter  closes  by  discussing  how  the 
performance  metrics  are  calculated  and  why  they  were  chosen. 

Chapter  4  provides  the  accuracy  and  confidence  data  produced  by  the  methodology. 
The  methodology  is  used  to  verify  the  functionality  of  the  NLO  algorithms  and  draw 
comparisons  between  their  performances  with  each  other  and  the  TA.  Chapter  5  begins  by 
describing  conclusions  drawn  from  the  results,  discusses  the  impact  of  these  conclusions, 
and  ends  by  describing  unexplored  questions  and  future  work. 
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II.  Background 


THe  mathematical  foundation  and  theory  underlying  this  thesis  is  presented  in  this 
chapter.  Section  2.1  defines  the  coordinate  systems  used  throughout  this  thesis  and 
describes  the  theory  behind  the  model  for  the  data.  Section  2.2  describes  the  mechanics 
behind  the  flightpath  of  the  Space  Vehicle  (SV).  The  Triangulation  Algorithm  (TA)  which 
is  used  for  an  initial  estimate  of  the  emitter’s  location  is  given  in  Section  2.3.  The  linear 
algebra  that  the  NLO  is  based  on  is  detailed  in  Section  2.4.  Section  2.5  presents  the  theory 
behind  the  technique  for  visualizing  the  confidence  in  the  NLO.  Monte  Carlo  Simulations 
(MCS)  as  a  means  of  visualizing  the  confidence  in  the  TA  are  given  in  Section  2.6.  This 
section  also  explains  how  MCS  are  used  as  truth  for  geolocation  confidence. 

2.1  Coordinates  and  Data 

The  standard  coordinate  system  that  is  used  throughout  this  research  is  Earth-Centered 
Earth-Fixed  (ECEF)  coordinates.  The  definition  of  ECEF  coordinates  is  described  [33]  as 
follows.  It  is  a  right-handed  coordinate  system  with  the  origin  at  the  center  of  mass  of  the 
Earth,  the  .r-axis  protrudes  at  the  zero  meridian,  the  z-axis  lies  along  the  earth’s  axis  and 
points  toward  the  North  pole,  and  the  y-axis  is  along  the  cross  product  of  the  z-axis  and 
.r-axis.  The  ECEF  coordinate  system  is  illustrated  in  Figure  2.1. 

The  LOS  data  is  described  by  two  spherical  angles  in  the  ECEF  coordinate  system. 
The  definition  of  these  angles  relative  to  the  ECEF  coordinate  system  is  shown  in 
Figure  2.2.  In  this  figure,  the  magenta  line  ij/  represents  an  arbitrary  Cartesian  LOS  vector 
as  given  by 
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The  yellow  line  represents  its  projection  onto  the  xy  plane.  Therefore,  6  is  the  azimuthal 
angle  measured  clockwise  about  the  z-axis  between  the  x-axis  and  the  projection  of  ij/  onto 
the  xy  plane.  4>  is  the  zenith  angle  measured  between  the  z-axis  and  i/r.  The  conversion  from 
spherical  to  Cartesian  coordinates  is  given  by 

if/x  =  sin(^)  cos  (8) 

ipy  =  sin(0)  sin(@)  (2.2) 

<AZ  =  co  s(0) 

Note  that  \\ifr\\  =  1  because  LOS  measurements  describe  only  the  direction  to  the  subject. 
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Figure  2.2:  Spherical  and  Cartesian  coordinate  convention 


The  conversion  from  ECEF  Cartesian  coordinates  to  spherical  coordinates  [21]  is 
given  by 

®  =  “''(£)  (23) 

^ 

cp  =  tan  - 

V 

It  is  arbitrarily  chosen  that  AOA  or  LOS  to  the  subject  is  measured  in  spherical  coordinates 
by  the  satellites.  If  the  positions  of  the  sensor  s  and  the  subject  x  are  given  as  Cartesian 
vectors  referenced  to  the  origin  of  the  ECEF  coordinate  system,  then  the  LOS  vector  is 
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simply  their  difference  vector.  Thus, 


<A  =  (2-4) 

||x  -  s|| 

There  are  imperfections  in  any  measurements.  Therefore,  the  AOA  measurements  are 
taken  to  be  random.  As  per  the  discussion  on  research  assumptions  in  Section  1.2.1,  a 
Gaussian  (or  Normal)  distribution  is  chosen  for  the  measurement’s  distribution.  Therefore, 
the  noisy  AOA  measurements  are  modeled  as 


(2-5) 

<P  ~  N  (filh  oj) 

In  these  equations  /i  represents  the  mean  and  cr  represents  the  standard  deviation  of  the 
respective  distribution.  It  is  assumed  that  there  is  no  bias  in  the  AOA.  Thus,  if  the 
subject  were  stationary,  He  and  fi(b  are  the  6  and  (b  angles  associated  with  if/  as  given 
in  Equation  (2.3).  With  the  coordinate  system  defined,  the  flightpath  of  an  SV  may  be 
described.  The  SV  described  in  the  following  section  is  used  to  evaluate  the  performance 
of  the  geolocation  algorithms. 

2.2  Space  Vehicle  Path 

To  use  a  geolocation  algorithm,  there  must  be  a  subject.  There  are  any  number  of 
subjects  that  might  be  used.  Thus,  a  particular  subject  was  chosen  for  simplicity.  The 
Space  Vehicle  (SV)  is  used  to  evaluate  the  performance  of  the  geolocation  algorithms.  As 
discussed  in  Section  1.2.1,  an  SV  was  chosen  as  the  subject  because  it  is  expected  to  be  the 
most  challenging  subject  for  geolocation. 
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An  SV  could  be  a  dead  satellite  or  space  junk.  These  may  be  in  orbit.  For  most 
orbits,  acceleration  is  nearly  constant.  To  make  this  subject  more  challenging,  a  reentry 
phase  was  added  to  introduce  rapid  deceleration  into  the  SV  model.  Realistically,  an  SV 
would  decelerate  due  to  friction  in  the  atmosphere  and  quickly  achieve  terminal  velocity. 
To  maintain  the  effects  of  acceleration  for  longer  periods  of  time,  the  SV  decelerates  until  it 
is  at  rest  at  its  crash  location.  Thus,  there  are  two  phases  to  the  SV’s  flightpath:  the  orbital 
and  reentry  phases. 

For  modeling  simplicity,  the  SV  flightpath  is  modeled  in  reverse.  The  SV  is  modeled 
in  two  phases.  There  is  an  acceleration  phase  and  an  orbital  phase.  These  are  modeled 
separately  and  put  together  at  the  end  of  this  section.  The  SV  flightpath  begins  with  its 
crash  location  and  then  accelerates  toward  its  orbital  phase.  The  orbiting  phase  is  modeled 
as  a  section  of  an  elliptical  path.  The  trajectory  is  simulated  using  the  mathematical 
development  that  follows.  Figure  2.3  illustrates  the  SV  flightpath  geometry.  The  flightpath 
is  simulated  in  the  xz- plane,  and  then  randomly  rotated  in  3-D. 

The  standard  form  of  an  ellipse  is  given  by 


This  equation  will  be  used  to  reach  the  parametric  form  of  an  ellipse  which  is  used 
throughout  the  rest  of  this  research. 

In  Equation  (2.6),  a  and  b  define  lengths  of  the  ellipse’s  axes,  the  largest  of  which 
defines  the  length  of  the  semi-major  axis.  The  smaller  value  is  the  length  of  the  semi-minor 
axis.  To  reach  the  parametric  form,  the  standard  form  of  an  ellipse  is  parameterized  as 

z  =  a  sin(0)  +  cz  (2.7) 

x  =  —b  cos  (0)  +  cx 
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e((j))  |  <h) 


x-axis 


Figure  2.3:  Geometry  for  space  vehicle  flightpath.  Note  that  for  visualization  purposes, 
this  figure  drastically  exaggerates  the  size  of  the  ellipse. 


Substituting  the  parameterized  variables  from  Equation  (2.7)  into  Equation  (2.6) 
produces 

(  a  sin(0)  +  cz  -  cz  (  -b  cos(0)  +  cx  -  cx 

\  «  )  +l  b 

sin2  cp  +  cos2  (f)  =  1 

This  equation  is  the  trigonometric  identity.  Therefore,  when  z  and  x  are  parameterized 
as  given  in  Equation  (2.7),  Equation  (2.6)  becomes  a  tautology.  Thus,  z  and  x  satisfy 
Equation  (2.6)  regardless  of  the  value  of  (p.  Therefore,  (p  is  a  parametric  variable  for  z  and 
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x.  An  ellipse  may  then  be  described  as  the  vector- valued  function 


e(0)  = 


—b  cos  <p  +  cx  a  sin  cp  +  cz 


(2.8) 


From  Figure  2.3,  the  center  of  the  ellipse  Ce  lies  along  the  z-axis.  Therefore,  cx  =  0 
and  Equation  (2.8)  is  simplified  to 


e(0)  = 


-b  cos  (p  a  sin  cp  +  Ce 


(2.9) 


This  elliptical  model  will  be  used  in  Chapter  3  to  simulate  the  SV  flightpath. 
Therefore,  the  ellipse’s  parameters  must  be  simulated.  The  following  discussion  describes 
how  these  parameters  are  established  in  the  simulations.  The  following  values  are  randomly 
generated:  Ce,  the  SV’s  apogee,  and  y.  Therefore,  one  of  the  ellipse’s  principle  axes  may 
be  found  via 


a  =  apogee  +  r  -  Ce 


(2.10) 


where  r  is  the  approximate  radius  of  Earth  and  is  given  the  value  6,378,100  meters. 

Next,  b  is  calculated.  However,  to  find  b.  [5  must  be  found  first.  The  z-componcnt  of 
e(0)|_ys  shown  in  Figure  2.3  is  used  to  find  yS  in  terms  of  other  known  quantities  as  given  by 


a  sin  -p  +  Ce  =  r  cos 


r 


y 

-a  sin/3  =  rcos  -  -  Ce 


Ce  -  r  cos  l 

sinyS  =  — - 2 

a 

Ce  -  r  cos  \ 

[3  =  arc  sin - - 

a 


(2.11) 
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j3  has  been  determined  in  terms  of  known  quantities.  Therefore,  the  value  of  b  may  be 
determined  by  using  the  ^-component  of  e(0)|_£  as  given  by 

-b  cos  -fi  = 
b 

Equations  (2.10)  and  (2.12)  provide  the  equations  needed  to  describe  the  ellipse  given  in 
Equation  (2.9).  Rather  than  parameterizing  the  flightpath  as  a  function  of  0,  the  flightpath 
must  be  described  as  a  function  of  time  t.  The  SV  impact  time  is  arbitrarily  defined  to  be 
at  t  -  0.  (Recall  that  the  SV’s  flightpath  is  simulated  in  reverse.)  The  point  of  impact  is 
chosen  to  be  at  e(0)|_^.  The  flightpath  is  then  modeled  as 

tT 

e(0  =  -b  cos  (ait  -  J3)  a  sin  (cot  -  /3)  +  Ce  (2.13) 

where  to  is  the  angular  velocity  of  the  SV.  The  flightpath  description  is  nearly  complete; 
however,  it  does  not  include  the  reentry  phase.  During  the  reentry  phase,  the  SV  starts  at 
rest  and  accelerates  until  it  reaches  its  orbiting  velocity.  Note  that  because  the  flightpath  is 
elliptical,  the  orbiting  velocity  slightly  oscillates.  The  altitude  of  a  typical  SV  at  the  end  of 
its  reentry  phase  is  approximately  known.  Therefore,  a  is  the  angle  that  satisfies 

cif  =  |e(cr)|  -  r  (2.14) 

This  results  in  a  quadratic  equation,  the  solution  to  which  is  given  in  the  code  given  in 
Appendix  A.  The  value  of  a  will  be  used  for  the  synthesis  of  the  reentry  phase.  The  reentry 
phase  is  given  in  mathematical  form,  and  will  be  followed  by  a  discussion  of  the  properties 
that  make  this  mathematical  model  attractive. 
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The  flightpath  during  the  reentry  phase  e,.(7)  is 


er(t)  = 


-b  cos  -  Pj  asin^j-cot  - p}  +  Ce 


A  a  +  fi 
where  tf  — 


co 


(2.15) 

(2.16) 


tf  is  defined  this  way  because  t  =  tf  when  cot  -  p  =  a.  A  realistic  orbital  velocity  v/  is 
more  important  for  realism  than  the  value  of  tf.  Therefore,  tf  is  left  as  a  variable,  and  a 
relationship  must  be  found  between  v/  and  tf.  The  value  of  tf  must  be  found  as  a  function 
of  Vf.  The  relationship  between  these  values  is  found  by  finding  the  velocity  of  the  SV  as  a 
function  of  time  vr(t),  ie 


vr(t)  =  |  —  er(f)| 
dt 


t  \2  (  t  x2 

\  2b co— sin  (cot  -  p)\  +  \2aco— cos  (cot  -  P) 

ff  /  \  ff 

-2co—  Jb2  sin2  (cot  -  p)  +  a2  cos1  (cot-  P). 

O' 


Vf  =  vb(tf)  =  2  co—  V^sin^a^-^cos^ 

O' 


a 


co 


vf 


2f  vfb2sin^a~Va2cos^ 


(2.17) 


This  equation  is  used  as  follows.  First,  a  is  found  by  solving  Equation  (2.14).  Next, 
Equation  (2.16)  is  used  to  linearly  relate  tf  and  co.  Finally,  co  is  found  in  terms  of  V/  by 
Equation  (2.17). 

The  orbital  phase  is  the  last  portion  of  the  flightpath  that  needs  to  be  defined.  Its 
definition  is 


eo(0 


-b  cos  (2 cot  -  cotf^j  a  sin  (2 cot  -  cotf )  +  Ce 


(2.18) 


This  definition  has  two  key  properties.  It  ensures  that  the  velocity  and  position  of  the  SV 
at  t  -  tf  are  the  same  for  both  phases  of  the  flightpath.  This  ensures  that  the  transition  is 
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continuous  for  position  and  velocity;  however,  higher  order  moments  such  as  acceleration 
and  jerk  will  not  be  constant.  Therefore,  the  final  equation  of  motion  for  the  S  V  flightpath 
is  given  by 


e(f)  = 


-b  cos  (jjtot  -  fi}  asva(j-<jot +  CeJ 
—b  cos  ^2  cot  -  totf^j  a  sin  (2  cot  -  cot  A  +  Ce 


t<tf 

t>tf 


(2.19) 


This  final  flightpath  equation  is  used  to  simulate  the  motion  of  the  subject.  The  next 
section  describes  a  typical  AOA  triangulation  method  for  geolocation. 


2.3  Triangulation  Algorithm 

If  extended,  LOS  measurements  from  several  satellites  should  nearly  intersect.  The 
Triangulation  Algorithm  (TA)  estimates  the  subject’s  location  as  the  point  where  these  LOS 
nearly  intersect  as  illustrated  in  Figure  2.4.  In  this  research,  the  TA  serves  two  purposes. 
First,  its  performance  is  compared  against  the  performance  of  the  NLO.  This  research 
proposes  the  NLO  as  a  means  of  mitigating  the  factors  that  the  TA  doesn’t  account  for  such 
as  system  geometry,  subject  motion,  and  time  delay.  Therefore,  the  NLO  should  provide 
improved  accuracy  and  confidence  in  its  geolocation  estimates  than  the  TA.  Second,  the 
NLO  requires  an  initial  guess  of  the  subject’s  position  to  begin.  The  TA  provides  the  initial 
guess. 

The  mathematical  foundation  for  the  TA  is  similar  to  that  given  in  [6]  and  [15].  The 
following  discussion  of  the  TA  has  been  presented  in  [5].  The  first  step  in  this  algorithm  is 
to  describe  the  minimum  distance  dj  to  the  zth  LOS  as  a  function  of  the  emitter’s  position  x. 
A  total  distance  D(x)  is  then  defined  as  given  by 

N 

D(x)  =  ^  d;2(x,  (2.20) 

1=1 
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Figure  2.4:  Triangulation  Algorithm  Notion.  The  Triangulation  Algorithm’s  geolocation 
estimate  minimizes  the  distances  to  the  lines-of-sight. 


where  w,  is  a  weighting  term  defined  on  [0, 1]  meant  to  describe  the  confidence  in  the  zth 
sensor’s  LOS  measurement. 

The  TA’s  geolocation  estimate  is  defined  as  the  point  in  space  that  minimizes  D. 
dj(x,  if/i)  is  quadratic  in  each  component  of  x.  D  is  also  positive  since  it  is  the  sum  of 
squared  values.  Therefore,  D  has  one  minimum.  This  minimum  is  the  value  of  x  that  is  the 
solution  to 

VD(x)  =  0.  (2.21) 

This  equation  results  in  a  system  of  three  linear  equations  in  three  unknowns.  This 
system  may  then  be  solved  by  any  standard  linear  algebra  technique  to  find  the  TA’s 
geolocation  estimate.  The  next  section  describes  the  mathematical  theory  underlying  the 
NLO. 
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2.4  NLO  Mechanics 


This  section  introduces  the  mathematical  theory  for  the  NLO.  NLO  is  a  means 
of  solving  systems  of  nonlinear  equations.  It  is  also  known  as  nonlinear  regression  in 
statistics.  Specifically,  the  Gauss-Newton  method  is  utilized  in  this  research  though  other 
approaches  could  conceivably  be  used.  NLO  begins  with  an  initial  approximation  for  the 
solution.  The  NLO  is  based  on  the  assumption  that  a  system  of  nonlinear  equations  may 
be  approximated  by  a  linear  system  in  the  neighborhood  of  the  solution.  The  linear  system 
is  solved  for  a  new  estimated  solution.  It  is  anticipated  that  the  estimated  solution  will  be 
closer  to  the  true  solution  to  the  nonlinear  system  of  equations  than  the  initial  guess. 

This  process  iterates  resulting  in  an  estimated  solution  that  converges  to  the  true 
solution  to  the  nonlinear  system.  In  this  way,  the  analytically  intractable  problem  of  solving 
the  system  of  nonlinear  equations  is  solved  to  a  desired  level  of  accuracy  by  iteratively 
solving  linear  systems  of  equations. 

The  following  discussion  will  introduce  NLO  with  a  simple  case  that  provides  valuable 
intuition  which  carries  forward  into  much  more  complicated  versions  of  the  problem.  The 
simple  case  is  a  classic  nonlinear  solver  known  as  Newton’s  method.  Newton’s  method  is 
developed  and  discussed  in  the  following  section  similar  to  the  explanation  given  in  [21]. 

2.4.1  Newton’s  Method  and  Intuition. 

Newton’s  method  was  originally  posed  as  a  means  of  finding  the  roots  of  polynomials 
of  any  order.  This  is  useful  because  there  is  no  analytical  means  of  finding  the  roots  of  a 
fifth  or  higher  order  polynomial.  Newton’s  method  is  applicable  to  more  than  polynomials. 
It  can  be  applied  to  most  differentiable  functions.  Therefore,  for  a  function  /(•),  Newton’s 
method  may  be  used  to  find  values  of  x  such  that  f(x)  =  0.  While  Newton’s  method  was 
originally  posed  for  finding  the  roots  of  functions,  it  is  suited  to  find  the  values  of  x  that 
produces  any  output  c  from  the  function.  Simply,  define  some  new  function  g(.v)  =  f(x)  -  c 
and  use  Newton’s  method  to  find  the  roots  of  g(-). 
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The  first  steps  in  Newton’s  method  are  illustrated  in  Figure  2.5.  The  objective  of 
Newton’s  method  is  to  find  the  input  xd  to  some  function  /(•)  that  produces  a  desired  output 
f{xd).  Therefore,  Newton’s  method  finds  the  value  of  x  given  by 

\x\f{x)  =  f(xd))  (2.22) 

There  may  be  several  values  of  x  for  which  /(•)  produces  f(xd ).  The  particular  solution 
that  is  sought  after  (the  desired  value  of  x)  is  xd. 


Figure  2.5:  Demonstrating  the  first  two  steps  in  Newton’s  method 


In  Newton’s  method,  both  f(xd)  and  an  initial  guess  at  xd  called  ,r0  are  known.  It 
is  assumed  that  /  is  approximately  linear  in  the  neighborhood  about  jc0.  Therefore,  the 
function  over  this  region  may  be  approximated  by 

f(x)  *  f'(x 0)(x  -  x0)  +  f(x 0)  (2.23) 

The  initial  guess  is  represented  by  the  blue  point  in  the  figure.  Next,  the  difference 
between  f(x)  and  f(xd),  typically  referred  to  as  the  residual  A/,  is  found.  Using  the  linear 
approximation  and  solving  this  equation  for  *  where  x  =  xd  should  produce  an  estimate 
of  xd  called  xd.  Equation  (2.23)  is  translated  into  point-slope  form  and  a  more  intuitive 
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notation,  ie 


(2.24) 


A,  df0 
A /  =  —  (*</  -  .x:0) 
ax 


Solving  this  equation  for  xd  yields 


a  fdx 
Xd  =  .To  +  A/  — 

dfo 


(2.25) 


f(xd)  should  be  closer  to  the  objective  value  f(xd)  than  f(x{)).  xd  is  represented  in  Figure  2.5 
by  the  yellow  point.  Next  let,  x\  =  xd.  This  process  now  repeats  by  letting  x\  be  the  initial 
guess  for  the  new  iteration.  This  is  illustrated  by  the  purple  line.  The  value  of  xd  will 
quickly  approach  the  solution  Equation  (2.22)  if  locally  linear.  This  iteration  stops  once 
the  residual  A /  has  reached  some  small  value  chosen  by  the  user. 

It  should  be  noted  that  this  technique  may  fail  as  a  result  of  using  an  inaccurate  starting 
point.  Two  cases  are  demonstrated  in  Figure  2.6. 


Figure  2.6:  Cases  where  Newton’s  method  fails.  The  top  starting  point  (the  triangle)  is  at  a 
peak,  so  the  slope  is  0  at  this  point.  The  lower  starting  point  (the  square)  will  converge  to 
the  black  point  rather  than  the  green  point. 


In  Figure  2.6,  the  yellow  triangle  and  square  are  initial  guesses,  the  green  point  is  the 
desired  point,  and  the  black  point  has  the  same  value  of  /(•)  as  the  green  point.  The  guess 
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at  the  local  maximum  in  this  plot  is  an  example  of  plateauing.  Here  the  slope  is  almost 
zero.  Thus,  the  jj  term  in  Equation  (2.25)  is  undefined  or  nearly  infinite. 

The  lower  initial  guess  is  another  example  of  a  inaccurate  starting  point.  Here 
Newton’s  method  will  converge  to  the  black  point  rather  than  the  green  point.  If  x,i 
represents  a  physical  value,  then  the  value  of  x  at  the  black  point  may  be  absurd.  For 
instance,  in  geolocation,  the  purple  point  may  represent  a  subject  with  a  speed  greater  than 
the  speed  of  light.  Convergence  to  an  undesired  solution  is  the  more  common  cause  of 
failure  in  Newton’s  method.  For  this  reason,  the  initial  guess  at  Xd  must  be  close  to  the  true 
value. 

All  of  the  properties  and  problems  regarding  Newton’s  method  have  analogous 
counterparts  in  NFO.  Therefore,  the  intuition  provided  by  Newton’s  method  carries 
forward  into  the  NFO  problem  where  systems  of  nonlinear  equations  are  considered.  Fike 
Newton’s  method,  NFO  approximates  a  nonlinear  system  of  equations  with  a  linear  system 
of  equations.  Finear  systems  of  equations  are  iteratively  solved  to  approximate  the  solution 
to  the  nonlinear  problem. 

2.4.2  Linear  Least  Squares. 

NFO  operates  on  an  overdetermined  system  of  nonlinear  equations.  Since  a  system 
of  nonlinear  equations  can  often  not  be  solved  analytically,  it  is  locally  approximated  by  a 
system  of  linear  equations  [30].  In  many  applications,  the  system  of  nonlinear  equations  is 
overdetermined.  That  is,  there  are  more  equations  than  unknowns. 

Therefore,  a  least  squares  solution  to  the  linear  system  is  used.  The  objective  of 
linear  least  squares  is  to  solve  the  following  system  for  x  [31].  Note  that  the  following 
development  of  the  linear  (and  weighted)  least  squares  solution  follows  the  development 
given  in  [31].  The  development  starts  with 


Ax  =  Q, 


(2.26) 
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where  A  is  m  x  n  where  m  >  n,  x  is  an  n  x  1  vector,  and  i>  is  a  m  x  1  vector.  In  the 
general  case,  where  m  may  be  greater  than  n  with  rank(A)  =  n  this  system  is  unsolvable. 
Therefore,  the  “closest”  solution  is  found.  The  closest  solution  is  defined  as  the  solution 
which  minimizes 

£2  =  ||n-Ax||2  (2.27) 

where  E  is  a  value  representing  the  error  in  a  solution  to  Equation  (2.26). 

The  value  of  x  that  minimizes  E2  is  the  least  squares  solution  to  Equation  (2.26).  From 
[31],  the  least  squares  solution  x  is  the  solution  to 

ATAx  =  ATQ.  (2.28) 

Note  that  ATA  is  invertible  if  and  only  if  the  columns  of  A  are  linearly  independent. 
This  is  true  for  most  cases  involving  random  data,  especially  if  m  »  n.  Therefore,  the 
least  squares  solution  x  is 

x  =  (ATA)_1ATn  (2.29) 

Some  systems  of  equations  may  be  attenuated  by  random  noise.  In  this  case,  the 
equations  in  the  system  may  not  be  equally  reliable.  For  this  situation,  the  best  solution  to 
Equation  (2.26)  is  the  solution  that  is  most  likely  to  minimize 

E2  =  ||Wfi  -  WAx||2  (2.30) 

where  W  is  a  matrix  of  weights  for  each  equation.  This  solution  is  known  as  the  weighted 
least  squares  solution  xw.  It  is  the  solution  to 

WAxw  =  WQ  (2.31) 
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The  solution  to  this  equation  is  given  by 


(WA)TWAxw  =  (WA)TWft 
AWTWAxw  =  ATWTWft 

xw  =  (ATWTWA)_1ATWTWft.  (2.32) 

Notice  that  the  symmetric  matrix  WTW  appears  twice  in  Equation  (2.32).  This  matrix 
has  an  important  stochastic  interpretation.  Suppose  that  the  elements  of  ft  are  normally 
distributed.  In  this  case, 

2T1  =  WTW  (2.33) 

where  IE1  is  the  inverse  of  the  covariance  matrix  £  of  ft  [8],  [31].  The  covariance  matrix 
is  of  the  form 

=  o-jO-j  (2.34) 

where  cr?  is  the  variance  in  the  zth  equation  of  Equation  (2.31).  Note  that  some  of  the 
randomly  attenuated  equations  may  be  interrelated.  Random  attenuation  in  one  equation 
may  affect  another.  This  is  accounted  for  by  the  covariance  terms  of  £.  The  covariance 
terms  are  E(/,  i  ±  j. 

Therefore,  in  the  presence  of  a  normally  distributed  ft.  Equation  (2.32)  may  be 
reduced  to 

xw  =  (AT2r1A)~1AT2r1ft  (2.35) 

This  equation  gives  the  weighted  least  squares  solution  to  Equation  (2.31).  Given  this 

discussion,  the  nonlinear  weighted  least  squares  problem  that  the  NLO  solves  may  be 
developed. 
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2.4.3  Weighted  Nonlinear  Least  Squares. 

NLO  is  a  method  for  approximating  the  weighted  nonlinear  least  square  solution  to  a 
system  of  nonlinear  equations.  Let  £2(x)  be  an  m  x  1  vector  of  m  nonlinear  functions  where 


ft(x)  = 


/l(Xl) 


flfa) 


X  = 


Xi 


*2 


(2.36) 

(2.37) 


In  this  problem,  a  desired  output  £2(xrf)  is  known.  A  particular  solution  xd  is  desired  given 
by 

{x|Q(x)  =  £l(xd)}  (2.38) 

where  xd  is  the  input  vector  that  produces  the  desired  output  flix,/).  Compare 
Equation  (2.38)  to  Equation  (2.22).  These  are  very  similar  except  that  fi(x)  represents 
a  set  of  functions. 

Because  the  system  of  nonlinear  equations  may  be  overdetermined,  there  may  not  be 
a  solution  to  Equation  (2.38).  Even  if  the  system  were  exactly  determined,  an  analytical 
solution  may  be  intractable.  Therefore,  the  objective  is  to  find  the  weighted  least  squares 
solution.  However,  since  the  system  of  equations  is  nonlinear,  even  the  weighted  least 
squares  solution  cannot  be  analytically  calculated  as  in  Equation  (2.35).  The  Gauss-Newton 
method  (which  is  called  NLO  in  this  thesis  due  to  the  application)  may  be  employed  to 
resolve  this  difficulty  [2]. 

In  the  region  about  a  some  initial  guess  of  xd  called  x0,  it  is  assumed  that 


£1(x)  ~  ft(x0)  +  J(x)(x  -  x0) 


(2.39) 
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where  J(x)  is  a  matrix  known  as  the  Jacobian  [2].  The  Jacobian  is  defined  as 


(V/!(X))T 

dfi  (x) 
dx\ 

4/i(x) 

0x2 

4/i(x) 

dxm 

J(X)  = 

(V/2(x))t 

= 

d/2(x) 

dx\ 

dx2 

dxm 

(V/;„(x))T 

df,„{x) 

dx\ 

dfm(x) 

dx2 

dfm(x) 
dxm  _ 

(2.40) 


Next,  Equation  (2.39)  is  expressed  in  an  intuitive  form  similar  to  that  given  in 
Equation  (2.24),  ie 

AO,  =  J(x0)(xf/  -  x0)  (2.41) 


where  AQ(x)  is  the  residual.  To  make  the  equation’s  intuitive  nature  clear,  let  Ax  =  xd  -  x0, 
so  the  equation  is  expanded  as 


A/i 
A  h 

Afm 


dJ^Axx  +  ^Ax2  +  •  •  •  + 
dJ^Axx  +  dJg*Ax2  +  •  •  •  + 

4^50)  AXl  +  ^Ax2  +  •  •  •  + 


m*0)  Ay 

f)  y  L*A-m 

d/2(xo)  . 
f)  r  AAA/77 
UXm 


dfn(X()) 

dxm 


Ax 


(2.42) 


From  this  equation,  it  can  be  seen  that  the  Jacobian  is  analogous  to  the  derivative  of  f(x )  in 
Equation  (2.23). 

Similarly  to  Newton’s  method,  the  objective  is  to  solve  Equation  (2.41)  for  x,/.  J(x0) 
is  an  m  x  n  matrix  where  m  >  n.  Often,  this  system  is  overdetermined  and  unsolvable.  Also 
note  that  weighting  has  not  been  included.  Therefore,  the  weighted  linear  least  squares 
solution  xd  to  Equation  (2.41)  is  the  solution  to 


WJ(xo)(x,  -  xo)  =  WA£1 


(2.43) 
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The  solution  to  this  equation  for  xd  as  given  by  Equation  (2.35)  is 


xd  =  x0  +  (J(x0)tI:-1J(xo))-1J(xo)t2:-1A£1  (2.44) 

As  in  Newton’s  method,  the  value  of  xd  should  produce  a  value  of  £l(xd)  that  is  closer  to 
the  desired  output  fl(xd).  Therefore,  xd  becomes  the  guess  for  the  second  iteration  of  the 
method  by  letting  Xi  =  xd.  The  process  then  repeats.  After  multiple  iterations,  the  value  of 
xd  converges  to  the  weighted  least  squares  solution  of  Equation  (2.38).  The  process  halts 
once  the  residual  reaches  a  value  specified  by  the  user. 

Recall  from  Section  2.4.2  that  in  the  presence  of  Gaussian  noise,  £  is  the  covariance 
matrix  of  that  distribution.  Furthermore,  it  is  shown  in  [18]  that  xd  converges  to  the  local 
maximum  likelihood  estimate.  That  is,  xd  converges  to  the  most  probable  x  given  £1  and 
En-  The  remaining  task  is  finding  a  means  of  describing  the  confidence  in  xd. 

Fortunately  this  method  naturally  produces  a  means  of  computing  the  confidence  in 
xd.  The  next  section  describes  how  NLO  is  used  to  calculate  the  confidence  and  produce 
an  associated  visualization. 

2.5  Estimating  and  Visualizing  Confidence 

Sections  2.3  and  2.4  provide  the  basis  for  the  geolocation  algorithms  developed  in  this 
thesis;  however,  the  visualization  scheme  has  not  been  addressed.  The  objective  of  this 
section  is  to  develop  a  method  for  visualizing  the  confidence  in  an  estimated  parameter  x. 
This  will  be  accomplished  with  the  covariance  matrix  of  x  called  Lx. 

The  covariance  matrix  describes  the  variance  and  covariance  in  x.  The  covariance 
matrix  needs  to  be  calculable  and  put  into  a  visually  intuitive  form.  The  visualization 
scheme  will  be  given  first,  followed  by  an  explanation  of  how  Lx  is  calculated. 

A  method  for  creating  this  visualization  is  found  in  [3].  This  visualization  method 
begins  by  calculating  the  Cholesky  Decomposition  of  the  covariance  matrix.  It  decomposes 
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a  matrix  £x  as  [31] 


Ex  =  LL*  (2.45) 

where  L  is  a  lower  triangular  matrix.  Note  that  the  Cholesky  Decomposition  only  operates 
on  Hermitian  positive-definite  matrices  [3].  Fortunately,  all  covariance  matrices  are 
Hermitian  and  positive  semi-definite  [31].  In  this  way,  L  is  analogous  to  the  square -root  of 
Sx- 

Using  L,  it  is  possible  to  generate  samples  of  x  that  follow  a  distribution  with  a 
covariance  matrix  of  Ex  from  normally  distributed  samples  y  as  given  by 


Y~  JV(0,I) 

x  =  LY  (2.46) 

— >x~  JV(  0,£x) 

This  procedure  will  create  a  scatter  of  x  values.  However,  if  rather  than  pulling  Y 
values  from  a  normal  distribution,  they  are  samples  of  a  sphere  of  radius  cr,  then  x  will 
sample  the  contour  of  its  distribution  associated  with  cr.  The  contours  of  Lx  are  ellipsoids 
[2].  Therefore,  the  process  for  visualizing  these  ellipsoids  is  given  by 

Y  ~  Sphere(Radius  =  cr) 
x  =  LY 

— >  Plot  Surface(x) 

where  cr  is  the  standard  deviation  associated  with  the  ellipsoid. 

The  standard  deviation  relates  to  a  percent  probability  that  the  samples  x  should  exist 
inside  of  the  ellipsoid.  Thus  the  percent  probability  pcr  that  x  will  exist  inside  of  an  ellipsoid 
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generated  with  a  given  cr  may  be  found.  This  ellipse  is  what  is  referred  to  in  this  thesis 
as  the  error  surface.  It  should  be  noted  that  the  equation  relating  p(T  and  cr  depends  on 
the  dimension  of  Lx.  It  is  commonly  known  that  in  a  1-D  normal  distribution  that  3cr 
corresponds  to  a  99%  probability.  This  is  not  true  for  a  3-D  normal  distribution.  For  a  3-D 
normal  distribution,  3<r  corresponds  to  a  probability  of  97%. 

This  procedure  provides  a  method  for  creating  a  visualization  of  the  covariance  matrix; 
however,  it  does  not  provide  a  means  of  finding  the  covariance  matrix.  Fortunately,  the 
covariance  matrix  can  easily  be  produced  for  the  NLO.  The  following  development  closely 
follows  that  given  in  [3];  however,  it  has  been  tailored  to  this  application.  First  the  definition 
of  Ex  is  given  by 

=  E[(xd  -  xd)(xd  -  xd)T]  (2.47) 

where  \,i  is  the  value  of  x  that  solves  Equation  (2.38)  exactly.  Recall  that  the  reason  xd 
cannot  be  exactly  found  is  because  of  the  random  noise  in  the  system.  In  Equation  (2.47), 
xd  is  a  random  variable  representing  the  weighted  nonlinear  least  squares  solution  given  by 
Equation  (2.44).  Next,  let  the  noise  in  the  system  be  described  by  a  random  variable  v  with 
the  properties 


E(v)  =  0 
E(vvT)  = 


(2.48) 

(2.49) 


Using  this  description  of  the  noise  in  the  system  and  ignoring  the  shift  x0,  the  left  hand 
side  of  Equation  (2.43)  is  rewritten  as 

JWWJxrf  =  JTWTWJxrf  +  JTWTWv.  (2.50) 
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The  only  random  attenuation  in  the  system  is  due  to  v.  Therefore,  the  covariance 
matrix  is  =  WTW.  Also  keep  in  mind  in  this  development  that  xd  and  xd  are  random 
variables.  Equation  (2.50)  is  then  simplified  as 

JTSn1J(x,-x,)  =  rSn'v.  (2.51) 

This  equation  will  be  used  to  find  a  simple  means  of  computing  in  terms  of  known 
values.  Lji  is  found  by  manipulating  both  sides  of  Equation  (2.51),  ie 

E  [(jTl£ J(x,  -  x,))  (ri^'jfx,  -  x,))T]  =  E  f(jTLn1v)  (jT^  v)"] 

E  [(j^1  J(xd  -  xc/))  ((x,  -  xd)TJTl£j)]  =  E  [(jT^v)  (vTEa1  j)] 

(jT2n  J)  E  [(xd  -  xd)(xd  -  x,)T]  (jTEn1  j)  =  r^E  [vvT]  I^1  J 
(jT2n1  J)  E  [(xd  -  xd)(xd  -  x,)T]  (jTEn1  j)  = 

(jT^n  J)  =  I 

Ex  =  (jT2n  j)  ‘  (2.52) 

This  equation  provides  a  simple  means  of  calculating  £x  from  the  known  matrices  J 
and  Zo.  Next,  the  location  of  the  error  surface  in  space  must  be  determined. 

To  this  end,  it  will  be  demonstrated  that  xd  is  an  unbiased  estimate  of  xd.  By  taking 
the  expected  value  of  both  sides  of  Equation  (2.51)  we  have 

E  [jT£a  J(N/  -  xrf)j  =  E  [jTSnv] 

JT^1J(E[xrf]-xrf)  =  J%!£[v] 

J(E  [xrf]  -  xd)  =  0  from  Equation  (2.48) 

=  E  [xd]  (2.53) 
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This  final  equation  follows  since  it  is  assumed  that  J  1  LoJ  is  invertible  and  therefore 
has  no  null  space  [31].  This  demonstrates  that  the  distribution  of  xd  is  centered  on  xd  and 
thus  is  an  unbiased  estimate. 

So  far,  this  visualization  method  is  only  applicable  to  the  NLO  because  no  method  has 
been  provided  to  calculate  Lx  for  the  TA.  Moreover,  no  method  has  been  given  to  evaluate 
if  the  covariance  matrix  of  xd  accurately  estimates  the  confidence  in  xd.  These  issues  are 
addressed  in  the  next  section. 

2.6  Monte  Carlo  Simulations 

In  this  section,  the  theory  and  utility  of  Monte  Carlo  Simulations  (MCS)  are  discussed. 
MCS  are  investigated  here  as  a  means  of  providing  the  confidence  surface  for  the  TA  and 
the  NLO.  The  confidence  surface  for  the  NLO  produced  by  MCS  is  treated  as  truth  for 
comparison  against  the  confidence  surface  produced  via  Equation  (2.52).  The  outputs  of 
some  systems  operating  on  random  variables  may  not  be  analytically  determined,  but  MCS 
provides  a  means  of  statistical  inference  from  simulations  [12],  [28]. 

The  fundamental  idea  behind  MCS  is  to  determine  the  distribution  of  a  system’s  output 
from  many  simulated  inputs  produced  from  the  inputs’  distributions.  In  our  case,  for  a 
given  sensor  configuration,  LOS  measurements  are  simulated  randomly  according  to  their 
distributions.  Each  set  of  random  measurements  is  used  to  estimate  the  subject.  This 
produces  a  distribution  of  subject  states.  This  method  is  illustrated  in  Figure  2.7. 

The  process  begins  by  using  the  N  sensors’  measurements  as  the  means  of 
their  respective  distributions.  The  random  measurements  resulting  from  the  sensors’ 
distributions  are  used  to  produce  geolocation  estimate.  This  process  is  repeated  until 
M  estimates  are  produced.  The  geolocation  covariance  matrix  is  estimated  from  the  M 
geolocation  estimates. 


35 


Figure  2.7:  Monte  Carlo  Simulation  covariance  matrix  estimation. 


The  MCS  produces  a  scatter  of  geolocation  estimates.  An  example  distribution  of 
position  estimates  is  shown  in  Figure  2.8.  Since  the  sensors’  measurements  are  used  as  the 
means  of  their  distributions,  the  scatter  of  geolocation  estimates  produced  by  the  MCS  will 
be  centered  on  the  estimate  given  by  the  sensor  measurements.  Thus  the  distribution  will 
have  the  same  mean  as  the  NLO. 

The  primary  purpose  of  producing  MCS  is  to  provide  truth  for  the  confidence  surface 
provided  by  the  NLO.  To  accomplish  this  task,  the  distribution  of  the  scatter  is  described 
in  terms  of  a  covariance  matrix. 
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Figure  2.8:  Scatter  of  geolocation  position  estimates  from  MCS. 


Let  the  number  of  MCS  estimates  be  M.  The  zth  estimate  is  x,-.  The  sample  covariance 
matrix  may  be  calculated  as  [13] 

1  M 

£x  =  y  Yj  (x'  “  E  [x'])  (x<  -  E  M)T  •  (2.54) 

1=1 

For  a  large  value  of  M,  the  sample  covariance  matrix  will  converge  to  the  true 
covariance  matrix.  Therefore,  the  sample  covariance  matrix  is  an  estimate  of  the  true 
covariance  matrix.  Observe  that  in  Figure  2.7,  the  method  for  producing  geolocation 
estimates  from  measurements  is  not  given.  MCS  can  be  used  to  produce  an  estimate  of  the 
true  covariance  matrix  for  any  geolocation  algorithm.  Therefore,  MCS  provides  a  means 
of  estimating  the  covariance  matrix  for  the  TA  and  NLO. 

Note  that  the  value  of  M,  for  which  the  sample  covariance  matrix  converges 
sufficiently,  is  unknown.  One  method  of  choosing  M  is  to  select  a  value  that  from 
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experience  has  been  shown  to  converge.  Another  method  is  to  keep  increasing  the 
number  of  MCS  samples  until  the  sample  covariance  matrix  changes  by  a  sufficiently 
smaller  amount.  This  provides  a  method  for  visualizing  the  confidence  estimates  for  the 
geolocation  algorithms. 

MCS  is  the  only  means  identified  presently  for  accurately  producing  the  covariance 
matrix  for  the  TA.  It  also  provides  the  ability  to  check  the  accuracy  of  the  covariance  matrix 
produced  by  the  NLO. 

The  application  of  these  algorithms  on  simulated  data  is  described  in  the  following 
chapter.  The  methodology  addresses  how  the  TA  and  NLO  are  compared  against  each 
other  and  how  the  NLO  may  be  improved  to  take  into  account  more  sources  of  error.  These 
and  other  practical  concerns  regarding  the  development  of  this  geolocation  method  and 
confidence  estimation  and  visualization  scheme  are  described  in  the  following  chapter. 
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III.  Methodology 


THis  chapter  discusses  the  methodology  for  evaluating  the  efficacy  of  Non-Linear 
Optimization  (NLO)  and  determining  which  NLO(s)  provides  the  best  performance. 
In  Chapter  1,  the  TA  is  described  as  a  standard  means  of  estimating  the  subject’s  location 
from  Line-of-Sight  (LOS)  measurements.  The  literature  review  in  Chapter  2  led  to  interest 
in  NLO,  otherwise  known  as  the  Gauss-Newton  method.  The  objective  of  this  research 
is  to  compare  the  performance  of  the  NLO  with  the  TA,  and  the  benefits  provided  by 
improvements  to  the  NLO.  Furthermore,  it  is  desirable  for  system  operators  to  be  capable 
of  visualizing  the  confidence  in  geolocation  algorithms.  Therefore,  a  comparison  of 
confidence  estimates  is  another  research  objective. 

The  methodology  used  in  this  research  is  to  simulate  scenarios,  generate  measure¬ 
ments,  and  run  the  various  algorithms  on  this  simulated  data.  The  first  stage  in  this  method 
is  simulating  the  scenario.  The  subject  being  localized  is  used  to  parameterize  the  sensor 
generation.  Therefore,  the  flightpath  generation  is  described  first. 

3.1  System  Geometry  Generation 
3.1.1  Time  Generation. 

The  measurement  times  are  the  first  part  of  the  scenario  simulation.  The  times  are  used 
to  generate  the  subject’s  flightpath  and  the  sensor  positions.  Therefore,  timing  generation 
is  the  first  stage  in  generation  of  the  system’s  geometry  (that  is,  the  positions  of  the  subject 
and  the  observing  sensors). 

Sensors  may  produce  LOS  measurements  of  the  subject  at  different  times.  Moreover, 
it  may  be  desirable  to  use  sensor  measurements  that  are  close  to  the  time  at  which  the  user 
would  like  to  geolocate  the  subject. 


39 


To  this  end,  measurement  times  are  generated  in  a  pseudorandom  manner.  The 
measurement  time  for  each  sensor’s  initial  measurement  time  t0  is  chosen  according  to 
a  uniform  distribution.  The  parameters  of  the  uniform  distribution  were  arbitrarily  chosen 
as 


t0-  Unif  (1,9)  s  (3.1) 

Subsequent  measurement  times  4  are  generated  in  the  same  manner,  using  the 
previous  time  as  a  starting  point  as 


4  ~  4-i  +  Unif  (1,9)  s  (3.2) 

This  process  is  performed  for  each  sensor  until  4  becomes  larger  than  a  threshold  /max. 
Sensor  measurements  and  truth  data  (the  subject’s  true  positions  and/or  velocities)  are 
produced  for  these  times.  The  time  values  determine  when  the  subject’s  path  is  sampled. 
In  the  next  section,  the  synthesis  of  the  SV’s  flightpath  is  described. 

3.1.2  Space  Vehicle  Generation. 

Given  that  the  performance  of  geolocation  algorithms  depends  on  the  flightpath  of 
the  subject  under  observation,  this  research  is  simplified  by  limiting  the  subject  to  one 
type.  The  most  challenging  subject,  the  Space  Vehicle  (SV),  is  used  to  evaluate  the 
performance  of  the  NLOs  and  TA.  The  model  for  the  SV,  detailed  in  Section  2.2,  contains 
sinusoidal  elements.  Sinusoids  are  infinitely  differentiable,  so  the  flightpath  contains  an 
infinite  number  of  moments.  Furthermore,  the  tremendous  speed  of  an  SV  exacerbates  the 
sources  of  error  discussed  in  Chapter  1.  Therefore,  an  SV’s  flightpath  seems  to  encapsulate 
the  greatest  extent  of  the  sources  of  geolocation  error. 
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The  flightpath  is  pseudorandomly  generated  by  randomly  producing  a  few  of  the 
flightpath  parameters  from  Section  2.2.  These  values  allow  the  remaining  parameters  to 
be  calculated. 

The  orientation  of  the  SV’s  flightpath  is  determined  by  randomly  generating  a  unit- 
vector  n  as  given  by 


n  ~  TV  (0, 13) 


n 


(3.3) 

(3.4) 


This  vector  is  orthogonal  to  every  vector  contained  in  a  plane.  The  flightpath  must  be 
contained  in  a  plane  perpendicular  to  this  plane.  This  setup  is  shown  in  Figure  3.1.  There 
are  an  infinite  number  of  planes  orthogonal  to  the  plane  shown  in  blue  in  Figure  3.1. 
The  plane  containing  the  flightpath  is  determined  via  a  uniformly  random  rotation  of  the 
flightpath  about  n.  The  MATLAB®  function  used  to  simulate  the  satellites  is  given  in 
Appendix  A.  The  shape  of  the  flightpath  is  created  by  randomly  generating  four  of  its 
parameters. 

The  arc-angle  (y  in  Figure  2.3),  apogee,  and  final  altitude  a  f  and  speed  at  the  end  of 
the  reentry  phase  v/  are  randomly  generated  as 


y  ~  Uni f  (5.5  x  10 6/r,ir)  Radians  (3.5) 

apogee  ~  Unif  (900  x  103, 1.5  x  106)  m  (3.6) 

af  ~  Unif(\50  x  103,450  x  103)  m  (3.7) 

vf  ~  Unif  (5  x  103, 7  x  103)  m/s.  (3.8) 


Equation  (3.5)  is  an  estimate  of  the  typical  ranges  these  values  would  have  for 
a  Low-Earth  Orbit  (LEO)  satellite.  The  remaining  three  parameters  are  the  result  of 
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testing  which  values  produce  reasonable  results.  The  approximate  orbital  speed  of  a 
LEO  satellite  is  7.6  km/s.  Therefore,  the  velocity  at  the  transition  to  the  reentry  phase 
is  randomly  produced  as  given  by  Equation  (3.8).  Its  apogee  is  approximately  1,120  km. 
Therefore,  Equation  (3.6)  keeps  the  apogee  near  this  value.  Equation  (3.7)  and  the  others 
in  combination  are  used  to  ensure  that  the  flightpath  is  never  inside  the  Earth.  The  result  of 
this  SV  generation  effort  is  given  in  Chapter  4. 


Figure  3.1:  SV  vehicle  flightpath  generation.  An  SV  flightpath  is  shown  in  yellow.  The 
vector  n  that  determines  the  orientation  of  the  flightpath  is  shown  in  magenta.  The  plane 
orthogonal  to  the  plane  containing  the  flightpath  is  shown  in  blue. 


Once  the  flightpath  has  been  generated,  positions  and  velocities  of  the  SV  along  the 
flightpath  may  be  calculated.  Positions  and  velocities  are  calculated  at  each  time  that  a 
measurement  is  taken  by  a  sensor.  This  is  done  to  provide  simulated  truth  data.  In  this 
way,  the  estimated  SV  state  at  these  times  may  be  compared  to  its  estimated  state.  The 
positions  are  calculated  according  to  Equation  (2.19).  The  velocity  is  calculated  via  the 
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first  derivative  of  Equation  (2.19).  The  next  section  will  describe  how  the  SV’s  flightpath 
is  used  to  parameterize  the  satellite  generation. 

3.1.3  Satellite  Generation. 

Though  the  algorithms  in  question  may  be  implemented  on  data  from  any  type  of 
sensor  that  measures  LOS,  satellites  are  simulated  for  this  research.  Satellites  are  chosen 
because  satellites  measure  the  LOS  to  the  subject  across  the  largest  distances.  However, 
some  satellites  may  also  be  very  close  to  the  subject  because  various  SV’s  may  approach  the 
lower  satellite  orbits.  Thus,  out  of  all  sensor  platforms,  satellites  provide  the  most  diversity 
in  range  to  the  subject.  The  first  step  in  the  satellite  generation  is  the  orbit  generation.  The 
code  that  performs  the  satellite  generation  is  given  in  Appendix  A. 

Lour  typical  satellites  orbits  were  considered:  Low  Earth  Orbit  (LEO),  Medium  Earth 
Orbit  (MEO),  Geosynchronous  Earth  Orbit  (GEO),  and  High  Earth  Orbit  (HEO).  There 
are  many  satellites  in  various  versions  of  these  orbits.  Therefore,  a  general  literature  review 
was  performed  to  develop  approximate  models  of  these  orbits.  While  the  models  may  not 
be  precise,  the  most  important  feature  of  the  orbits  for  this  research  is  their  orbital  distance. 
The  models  for  these  orbits  provide  the  desired  diversity  in  the  oribtal  distances.  The  model 
of  each  satellite  is  given  in  the  orbital  generation  function  OrbGen.m  shown  in  Appendix 
A.  The  orientations  of  the  orbits  are  chosen  pseudorandomly. 

The  orbits  are  generated  at  random  with  the  constraint  that  all  satellites  must  have 
LOS  to  the  SV.  In  addition  to  realism,  this  constraint  was  added  to  prevent  satellites  from 
observing  a  subject  over  unrealistic  distances.  For  example,  if  a  LEO  satellite  is  on  the 
opposite  side  of  the  Earth  relative  to  the  observed  subject,  the  error  in  its  LOS  propagates 
over  a  larger  distance  than  would  otherwise  be  possible.  This  constraint  is  implemented 
by  taking  advantage  of  the  randomly  generated  vector  n  from  the  previous  section.  The 
constraint  is  illustrated  in  Figure  3.2. 
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Figure  3.2:  Satellite  orbit  synthesis.  As  before,  the  magenta  line  represents  n.  The 
restriction  on  the  satellite’s  orbits  is  that  they  must  be  on  the  same  side  of  the  blue  plane  as 
the  SV. 


A  simple  means  of  ensuring  that  each  satellite  has  LOS  to  the  subject  at  the  time 
of  each  measurement  was  not  found.  Moreover,  this  requirement  is  unnecessary.  The 
purpose  of  this  requirement  is  to  ensure  that  sensors  do  not  observe  the  subject  at  absurd 
distances.  Therefore,  this  requirement  was  relaxed  and  substituted  for  the  requirement  that 
the  satellites  and  the  subject  are  on  the  same  side  of  the  Earth.  This  process  is  tedious  and 
is  left  to  the  code  given  in  Appendix  A.  A  notional  explanation  of  the  satellite  generation 
is  described  below. 

The  satellite  orbits  are  randomly  generated.  The  length  of  time  needed  to  observe  the 
subject  is  used  to  calculate  a  section  of  the  satellite’s  orbit  over  which  the  satellite  observes 
the  subject.  This  section  is  then  rotated  (if  needed)  until  both  the  endpoints  of  the  section 
are  on  the  same  side  of  the  Earth  as  the  subject.  The  sensors  are  on  the  same  side  of  the 
Earth  meaning  that  the  satellites  are  on  the  same  side  of  the  plane  whose  span  is  orthogonal 
to  n.  This  is  the  plane  shown  in  blue  in  Figure  3.2.  The  satellite  orbits  and  SV  flightpath 
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form  the  basis  for  LOS  generation.  The  method  for  simulating  the  LOS  for  each  algorithm 
is  described  in  the  section  3.2. 

3.2  LOS  Generation 

This  section  discusses  how  LOS  measurements  are  simulated.  LOS  measurements  are 
generated  for  three  purposes.  The  first  is  to  evaluate  the  comparative  performance  of  the 
NLO  to  the  TA.  The  second  is  verify  the  theoretical  performance  of  the  algorithms.  Lastly, 
the  LOS  data  is  produced  to  determine  the  benefits  provided  from  adding  complexity  to  the 
NLO.  Each  NLO  takes  into  account  certain  sources  of  error.  The  LOS  synthesis  for  the 
Static  NLO  is  described  in  the  next  section. 

3.2.1  Stationary  Subject. 

Each  NLO  includes  certain  sources  of  LOS  error  in  their  model.  The  simplest  NLO, 
the  Static  NLO,  only  takes  into  consideration  the  geometry  of  the  system.  It  does  not 
include  sensor  synchronization,  the  subject’s  velocity,  and  time  delay  in  its  model.  Sensor 
synchronization  means  that  the  sensors  obtain  their  measurements  at  the  same  time.  The 
time  delay  is  the  time  it  takes  for  a  signal  coming  from  the  subject  to  be  received  by  the 
sensors. 

The  Static  NLO  assumes  that  the  subject  is  stationary  over  the  duration  of  the 
measurements  used  to  provide  an  estimate.  This  assumption  is  inaccurate  for  fast  moving 
subjects  and  unsynchronized  sensor  measurements.  In  such  a  case,  the  NLOs  confidence 
surface  is  expected  to  be  overly  confident. 

To  verify  that  the  static  NLO  performs  correctly,  a  system  is  created  where  the  sensors’ 
measurements  are  perfectly  synchronized,  and  the  subject  is  stationary.  The  Static  NLO  is 
optimized  for  this  type  of  scenario.  If  the  static  NLO  has  been  correctly  implemented,  then 
its  confidence  estimate  will  be  ideal  and  should  match  its  MCS. 

Lor  the  static  NLO  the  LOS  generation  is  straightforward.  The  measurement 
distribution  follows  from  Ligure  3.3.  Since  the  subject  is  stationary,  the  ideal  LOS 
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measurement  is  merely  the  difference  vector  between  the  subject  and  the  sensor.  Therefore, 
since  it  is  assumed  that  the  measurements  have  no  bias,  the  means  of  the  6  and  <p 
distributions  are  the  spherical  angles  associated  with  this  difference  vector.  The  difference 
vector  if/  =  x-  s.  Given  this  difference  vector  and  taking  note  of  Equation  (2.3),  this  means 
the  distributions  corresponding  to  the  difference  vector  between  the  subject  and  a  sensor 
can  be  calculated  as 

H„  =  arctan  2  (if/y,  if/xj  =  arctan  2  (xv  -  sv,  xv  -  sT)  (3.9) 

AG  =  arctan  2  (  ^ jif/2x  +  if/],,  if/z j  =  arctan  2  ^  ^(x.r  -  sA-)2  +  (xv  -  sv)2,  xz  -  S- j 

where  He  and  are  the  respective  means  of  the  6  and  <p  distributions.  The  distributions 
are  assumed  to  be  Gaussian  as  given  by  Equation  (2.5).  Because  Equation  (3.9)  allows  the 
means  to  be  calculated,  only  the  variances  are  unknown. 

The  variances  of  the  6  and  0  distributions  are  determined  by  the  confidence  in  the 
LOS-measuring  sensor.  The  variances  are  arbitrarily  selected  in  this  research.  In  the  next 
section,  this  problem  is  expanded  to  include  subject  motion. 

3.2.2  Adding  Object  Motion. 

In  this  section,  the  LOS  data  takes  into  account  linear  motion.  The  subject  has  constant 
velocity,  and  the  time  delay  is  not  included.  The  utility  of  this  data  is  in  verifying  the 
functionality  of  the  Velocity  NLO.  The  Velocity  NLO  assumes  that  the  subject  moves 
in  approximately  a  straight  line  over  the  duration  of  the  measurements  used  to  provide 
a  geolocation  estimate,  and  it  does  not  include  time  delay  in  its  model.  Therefore,  the 
Velocity  NLO  should  converge  to  the  ideal  estimate  of  the  subject’s  state  and  provide  an 
ideal  confidence  surface. 

Because  the  time  delay  is  not  included,  the  sensors  measure  the  LOS  to  the  subject 
instantaneously.  Therefore,  the  means  of  the  sensor  distributions  are  calculated  in  the  same 
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Figure  3.3:  Static  NLO  LOS  synthesis.  The  coordinates  to  the  sensor  are  given  by  the 
vector  s.  The  subject’s  coordinates  are  given  by  the  vector  x.  if/  is  the  difference  vector 
between  the  subject  and  the  sensor.  6  and  (p  are  the  spherical  angles  associated  with  ij/. 


manner  as  in  the  static  case  described  by  Equation  (3.9).  The  only  difference  in  this  case  is 
that  the  subject’s  position  changes  with  time.  The  subject’s  position  at  a  given  time  p(t)  is 
described  by  the  vector- valued  function  as  given  by 

P(0  =  P(Olf=o  +  PF  (3.10) 

In  this  equation,  the  subject’s  velocity  is  the  first  derivative  of  the  position  with  respect 
to  time  and  denoted  by  p.  Note  that  the  velocity  is  a  constant  and  thus  not  a  function  of 
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time.  Different  velocities  are  used  to  test  the  NLOs.  Next,  the  time  delay  is  added  to  the 
linear  model. 

3.2.3  Adding  Time  Delay. 

The  LOS  data  generation  scheme  discussed  in  this  section  is  used  to  verify  the 
functionality  of  the  most  robust  NLO:  the  time  delay  and  velocity  NLO.  The  LOS 
measurements  described  in  this  section  include  time  delay  for  subject  with  a  linear 
flightpath. 

The  time  delay  is  particularly  relevant  for  satellite  systems.  For  example,  a  HEO 
satellite  may  have  an  approximate  orbiting  radius  of  60  x  106  m.  Radiation  from  a  subject 
sitting  on  the  Earth’s  surface  requires  approximately  (60  x  106  -  6.731  x  106) / 3  x  108  s 
«  200  ms  to  be  received  by  a  HEO  sensor.  Although  small,  this  time  delay  may  be 
significant  for  fast-moving  subjects.  The  dead  LEO  satellite  in  Section  3.1.2  has  a 
maximum  speed  of  7.6  x  103  m/s.  In  200  ms,  such  a  satellite  could  traverse  1.5  km!  Thus, 
by  the  time  the  sensor  has  observed  the  subject,  it  has  moved  1 .5  km.  The  time  delay  means 
that  sensors  are  not  observing  the  subject  where  it  currently  is  but  where  it  used  to  be. 

Therefore,  the  means  of  the  6  and  (p  distributions  are  determined  by  the  subject’s 
position  at  the  time  that  the  sensor  measures  it  minus  the  time  delay.  To  solve  this  problem, 
the  following  definitions  are  made.  The  time  delay  is  give  the  variable  td.  The  time  at 
which  the  subject  is  observed  by  a  sensor,  the  measurement  time,  is  called  tm.  The  LOS 
for  a  given  sensor  is  then  determined  by  the  subject’s  position  at  tm  -  td.  This  is  concerned 
with  data  simulation.  Therefore,  the  subject’s  motion  is  known.  Its  position  at  tm  is  known, 
and  the  sensor’s  position  is  also  known;  therefore,  finding  the  subject’s  position  at  tm  -  td 
requires  knowledge  of  only  td.  A  method  for  calculating  td  is  developed  using  Figure  3.4. 

The  time  delay  is  calculated  by  finding  an  equation  for  the  time  delay  in  terms  of  the 
the  subject’s  position  and  velocity  and  the  sensors  position.  Because  this  is  the  synthesis 
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Figure  3.4:  LOS  with  time  delay  synthesis  setup.  The  sensors  position  at  the  measurement 
time  tm  is  given  by  s(tm).  Two  of  the  subject’s  positions  are  shown  as  the  asterisks.  Its 
flightpath  is  shown  as  the  straight  lines  through  the  asterisks.  The  distance  between  the 
sensor  at  time  tm  and  the  subject  at  time  tm  -  td  is  denoted  by  d. 


problem,  the  subject’s  motion  is  known.  Therefore, 


x(Li  td)  —  *(tln )  (3.11) 

Because  the  subject’s  position  and  velocity  are  known,  the  unknowns  in  this  equation 
are  td  and  \(tm  -  td).  Another  independent  equation  is  developed  to  take  advantage  of  the 
fact  that  the  sensor’s  location  is  known,  yielding  an  equation  for  td. 

The  time  delay  is  the  time  it  takes  for  the  signal  emitted  from  the  subject  at  tm  -  td  to 
reach  the  sensor  at  tm.  The  signal  traverses  a  distance  of  d  at  a  rate  of  the  speed  of  light  c. 
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Therefore,  from  Figure  3.4, 


d  |x(rm  td )  s  (tn 


(3.12) 


Because  the  speed  of  light  is  a  known  constant,  td  may  be  given  as 


td  - 


l|x(tm  td )  S  (ttl 


(3.13) 


Substituting  Equation  (3.11)  into  the  above  equation  and  multiplying  by  c  produces 


ctd  —  \x(lm)  xtd  s 


(3.14) 


The  unknown  in  this  equation  is  td.  Squaring  both  sides  of  the  above  equation  produces 
a  quadratic  equation  for  td.  This  equation  may  be  solved  for  td.  Because  the  equation  is 
quadratic,  two  solutions  exist;  however,  only  one  is  realistic.  By  observing  the  nature  of 
the  solution,  the  correct  equation  for  td  is  selected.  One  of  the  solutions  is  negative.  The 
negative  solution  produces  a  value  of  tm  -  td  which  is  greater  than  tm.  By  the  problem’s 
setup,  the  value  of  td  must  be  positive,  so  the  negative  value  cannot  be  the  solution. 

Given  td,  the  means  of  the  6  and  0  distributions  may  be  calculated  as 

Hq  —  arctan  2  (xy(tm  t d )  sY(t/n ) ,  xx(tfn  !  d )  sv(tfn )  j  (3.15) 

Us  —  arctan  2  (xz(tm  t d )  —  sz(tm),  -^(Xj(7m  tj)  sx(7m))“  +  (Xy(t m  td)  Sy(tm))“| . 


With  td  known,  x(tm  -  ld)  is  also  known  as  per  Equation  (3.11). 

The  next  section  discuses  the  synthesis  of  LOS  data  for  a  model  of  the  most 
challenging  subject.  The  most  challenging  scenario  for  geolocation  is  considered  to  be 
the  case  of  observing  a  fast-moving  subject  such  as  an  SV. 
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3.2.4  Space  Vehicle  Scenario. 

Sections  3.2.1  to  3.2.3  describe  LOS  measurement  synthesis  for  unrealistic  subjects 
that  are  stationary  or  have  linear  flightpaths.  This  section  describes  LOS  data  generation 
for  the  most  challenging  subject:  the  Space  Vehicle  (SV). 

The  time  delay  cannot  be  analytically  calculated  for  the  SV.  Therefore,  time  delay  is 
calculated  as  part  of  the  subject’s  flightpath.  This  process  is  illustrated  by  Figure  3.5.  The 
idea  is  to  first  generate  the  measurement  times  for  the  sensors  as  described  in  Section  3.1.1. 
Next,  the  SV’s  flightpath  is  generated.  The  time  values  and  the  SVs  flightpath  parameters 
are  used  to  generate  sensor  positions.  The  time  delay  that  would  exist  between  the  sensors 
and  the  subject  at  the  time  of  measurement  tm  is  then  calculated  via 


td 


V  t m  )  S(tm)|| 

C 


(3.16) 


The  time  delay  for  each  measurement  time  is  added  to  the  measurement  time.  The 
SV  flightpath  is  simulated  at  the  time  corresponding  to  the  sensors’  measurements.  These 
positions  are  the  locations  of  the  SV  at  measurement  times  tm.  With  the  subject’s  position 
at  the  measurement  times  calculated,  the  LOS  angles  associated  with  these  positions  may 
be  simulated.  This  concludes  discussion  on  the  synthesis  of  data  for  the  most  challenging 
scenario. 


3.3  Static  NLO 

This  section  describes  how  the  static  NLO  is  constructed.  For  the  static  NLO,  it 
is  assumed  that  the  subject  is  stationary  while  sensors  produce  the  measurements  used 
to  geolocate  it.  As  discussed  in  Section  2.4,  at  the  core  of  this  method  is  the  linearity 
assumption. 

This  assumption  is  captured  in  the  Jacobian  of  the  measurements  as  a  function  of  the 
subject’s  state  £2(x).  This  function  outputs  the  ideal  (noiseless)  measurements  that  are 
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Figure  3.5:  Geometry  for  space  vehicle  LOS  synthesis.  The  subject’s  flightpath  is  shown 
as  the  curved  line.  Its  positions  at  the  measurement  times  are  illustrated  by  the  points  along 
the  path.  The  time  delay  associated  with  distance  between  the  sensor  at  time  tm  and  the 
subject  at  the  time  of  measurement  tm  is  denoted  by  t,t-  The  squares  illustrate  the  SV’s 
location  at  the  tm. 


associated  with  the  input  subject  state.  Equation  (3.9)  may  be  used  to  find  this  function 
for  one  pair  of  6  and  (p  measurements  in  the  static  scenario  (the  subject  is  stationary).  This 
equation  gives  the  means  of  the  6  and  (p  distributions,  respectively.  The  means  of  these 
angles  are  the  angles  that  should  be  measured  in  the  absence  of  noise.  Therefore,  the  ideal 
angles  that  should  be  measured  for  the  static  scenario  may  be  calculated  as 

iff  =  \  —  s 

6  =  arctan  2  (3.17) 

4 P  =  ^tan  2  (  ^  +  f2r  t/r,j 

where  s  is  a  vector  representing  the  sensor’s  position. 

To  localize  the  subject,  several  measurements  must  be  used.  Therefore,  the 
Jacobian  will  be  the  concatenation  of  the  Jacobian  of  each  measurement  as  described  in 
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Equation  (2.39).  Thus  the  Jacobian  at  k  -  1  is  given  as 
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(3.18) 


where  each  angle’s  subscript  indicates  the  measurement  with  which  it  is  associated.  The 
partial  derivatives  from  which  the  Jacobian  is  constructed  are  the  partial  derivatives  of  0 
and  (p  as  given  in  Equation  (3.17).  The  values  of  these  derivatives  at  x^_i  are  not  equal 
for  different  measurements  because  the  sensor  location  is  different  for  each  measurement. 
These  derivatives  are  given  compactly  in  the  form  of  the  gradients  of  the  measurement 
angles,  ie 


ll^l  I2  =  (X.V  -  Sx)2  +  (Xy  -  S  yf  +  (X,  -  S,)2 
ll«M2  =  (x.v  -  Sx)2  +  (xy  -  sv)2 


vx0 


Vx0 


(V'-Sr)(X;-S;)  (xr-sr)(X;-S'-) 

II<AII2II<Aa-vII  jWPlIiM 


(3.19) 


The  Jacobian  described  above  may  then  be  input  into  Equation  (2.44)  to  produce  an 
updated  estimate  of  the  subject’s  state  as  given  by 


X,  =  x*_!  +  ( Jj_ ,  £  ~ 1  Jj- 1 )~ 1 ,  E,,1  AH 


(3.20) 
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This  process  may  be  iterated  until  it  converges.  The  process  is  said  to  converge  when  the 
||Xjt  -  Xfc-i||  <  1  x  1  (T6 .  The  residual  term  Ail  is  the  difference  between  the  measured  LOS 
and  the  LOS  that  should  have  been  ideally  measured  if  the  subject’s  state  is  x*_i.  This  is 
expressed  by 


Ail  =  Q-n(x*_,)  (3.21) 

where  il  is  the  measured  LOS  angles.  This  process  requires  an  initial  estimate  to  begin 
using  the  NLO  at  k  =  0.  The  initial  estimate  is  provided  via  triangulation  as  described  in 
Section  2.3.  The  remaining  unknown  quantity  from  Equation  (3.20)  is  the  matrix  I£i. 

I,q  is  the  covariance  matrix  of  the  measurements.  The  covariance  matrix  is  arbitrarily 
defined  as 


Sn  =  (T-l  (3.22) 

where  the  standard  deviation  cr  =  5  x  10-6  rad.  (In  Chapter  4  the  change  in  geolocation 
performance  resulting  from  difference  standard  deviations  is  discussed).  With  the 
covariance  matrix  defined,  all  of  the  terms  in  Equation  (3.20)  have  been  described. 

The  covariance  matrix  for  the  measurements  may  be  combined  with  the  Jacobian  via 
Equation  (2.52)  to  produce  the  covariance  matrix  of  the  subject’s  state  as  given  by 


Ex  =  (jT£n'j)"'  (3.23) 

If  the  initial  state  estimate  is  close  to  the  true  state,  this  NLO  will  converge  to  an 
unbiased  least-squares  solution  to  the  static  geolocation  problem.  The  next  section  will 
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address  the  construction  of  the  Velocity  NLO.  The  Velocity  NLO  optimizes  for  the 
subject’s  motion. 

3.4  Velocity  NLO 

In  this  section,  an  NLO  is  developed  to  solve  the  geolocation  problem  where  the 
subject’s  position  changes  significantly  over  a  short  amount  of  time.  If  all  of  the  sensors 
measure  the  LOS  to  the  subject  at  the  same  time  and  subject’s  state  is  estimated  at  that  time, 
then  the  subject’s  velocity  is  irrelevant.  In  this  case  the  geolocation  problem  simplifies  to 
the  static  scenario.  The  more  interesting  problem  is  the  scenario  where  the  sensors  do  not 
measure  LOS  at  the  same  time.  In  this  case,  a  time  at  which  to  estimate  the  subject’s  state 
must  be  selected. 

We  let  the  time  when  a  sensor  produces  a  LOS  measurement  of  the  subject  be  tm. 
The  time  at  which  the  user  desires  to  estimate  the  subject  is  te.  For  this  type  of  LOS  data, 
it  is  assumed  that  the  subject’s  motion  is  approximately  linear  over  short  spans  of  time. 
Therefore,  the  subject’s  motion  is  described  similar  to  Equation  (3.10)  by 

p(0  =  p(OUo  +  P  t  (3.24) 


where  the  subject’s  position  is  given  by  p(t)  and  its  velocity  is  given  by  p.  As  before,  the 
velocity  is  constant.  There  is  now  an  additional  unknown  in  this  problem:  the  velocity  of 
the  subject.  Therefore,  the  subject’s  state  now  includes  velocity  and  is  ordered  as  given  by 


x  = 


X,  X.v  Xy 


(3.25) 


As  before,  the  NLO  requires  that  LOS  measures  are  described  in  terms  of  the  subject’s 
state  by  a  function  herein  called  fi(x)  .  In  this  case,  the  NLO  is  being  used  to  estimate  the 
subject’s  state  at  an  estimation  time  which  may  or  may  not  be  the  same  as  the  measurement 
time.  This  relationship  is  developed  in  the  manner  illustrated  in  Ligure  3.6.  The  subject’s 
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position  is  projected  in  time  via  its  model  of  motion  to  measurement  time  from  the 
estimation  time.  The  relationship  between  the  measurements  and  the  subject’s  state  is  then 
identical  to  that  of  the  static  case. 


Figure  3.6:  Setup  for  the  velocity  NLO.  The  subject’s  velocity  is  used  to  project  its  state 
from  the  estimation  time  to  the  measurement  time. 


From  Figure  3.6,  it  can  be  seen  that  the  difference  between  the  subject’s  position  at 
the  measurement  time  and  its  position  at  the  estimation  time  may  be  calculated  as 


P (tm)  =  p (te)  -  (te  -  tm) p.  (3.26) 

Because  the  objective  is  to  estimate  the  subject’s  state  at  te,  the  true  LOS  measurements 
associated  with  the  subject’s  state  at  te  must  be  calculated.  If  calculable,  this  allows  the 
LOS  measurements  to  be  compared  against  the  true  LOS  that  should  be  measured  if  the 
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estimated  subject  state  were  its  true  state.  The  difference  between  the  measured  LOS  and 
the  true  LOS  associated  with  the  estimate  of  the  subject’s  state  is  the  residual  term  Afl. 

It  is  assumed  that  in  a  general  scenario,  the  estimation  time  does  not  coincide  with  a 
measurement.  Equation  (3.26)  gives  an  equation  for  the  subject’s  state  at  the  estimation 
time  in  terms  of  the  measurement  time.  Therefore,  the  ideal  measurements  are  a  function 
of  the  subject’s  state  at  the  measurement  time  in  terms  of  its  state  at  the  estimation  time. 
This  is  mathematically  given  by 


£1  (x(fm)|x(fe))  =  £2 


VL 


P(L )  (L  tm)p 

p 


v 


(3.27) 


As  in  the  static  case,  the  NLO’s  6  and  0  measurement  angles  depend  only  on  the 
subject’s  position  at  the  measurement  time.  However,  because  the  measurements  are  now 
described  in  terms  of  the  subject’s  state  (position  and  velocity  in  this  case)  at  the  estimation 
time,  the  subject’s  velocity  must  also  be  taken  into  account.  Therefore,  Equation  (3.17)  is 
modified  for  the  scenario  with  a  constant- velocity  subject,  ie 


iff  =  (p (te)  -  p (te  -  tm ))  -  S 

6  =  arctan  2  (if/y,  if/^  (3.28) 

0  =  arctan  2  (  ^jif/2x  +  iff2,  iftz  j 

More  than  one  LOS  measurement  is  used  for  geolocation,  so  Q  (x)  is  the  concatenation 
of  6  and  0  and  functions  of  the  subject’s  state  for  each  sensor.  This  concatenation  of 
functions  is  linearly  approximated  using  the  Jacobian.  As  before,  the  Jacobian  contains 
the  first  partial  derivatives  of  6  and  0  with  respect  to  each  variable  in  the  subject’s  state  for 
each  sensor  at  the  current  estimate  of  the  subject’s  state.  In  Equation  (3.18)  the  Jacobian 
contained  only  partial  derivatives  with  respect  to  the  subject’s  position.  Now  the  subject’s 
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state  includes  its  velocity,  and  so  the  Jacobian  also  contains  the  partial  derivatives  of  6  and 
(p  with  respect  to  velocity.  The  Jacobian  for  the  velocity  NLO  is  given  by 
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These  partial  derivatives  are  taken  with  respect  to  6  and  (p  as  given  in  Equation  (3.28). 
The  result  of  this  differentiation  is  similar  to  that  given  by  Equation  (3.19);  however,  the 
derivatives  are  too  lengthy  to  be  written  in  this  document.  The  derivatives  are  recorded  in 
VNLO_NTD.m  given  in  Appendix  A. 

The  covariance  matrix  of  the  measurements  Zo  is  defined  just  like  in  the  static 
case.  The  NLO  is  then  used  to  improve  estimates  of  the  subject’s  state  according  to 
Equation  (3.20). 

The  NLO  requires  an  initial  estimate  of  the  subject’s  state.  Previously,  the  initial 
estimate  contained  only  position.  Now  the  initial  estimate  includes  velocity.  The  initial 
guess  of  velocity  is  produced  by  interpolating  between  initial  position  estimates.  The 
position  estimates  are  provided  via  triangulation.  The  initial  velocity  estimate  is  calculated 
by 


Pofe,)  = 


Po(tei)  ~  Mtei-t) 


tP  ~  t, 


<?/- 1 


i  >  1  (3.30) 


where  tgj  is  the  zth  estimation  time  and  the  0  subscript  represents  that  this  is  the  initial 
velocity  estimate  which  is  used  to  start  the  NLO.  Note  that  because  there  is  no  estimation 
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prior  to  the  first  estimation  time,  the  velocity  at  the  first  estimation  time  cannot  be 
interpolated.  Therefore,  it  is  assumed  that  the  velocity  at  the  initial  estimation  time  is 
the  velocity  at  the  subsequent  estimation  time. 

Since  the  NLO  also  optimizes  the  subject’s  velocity,  the  covariance  matrix  of  the 
subject’s  velocity  is  also  produced.  The  covariance  matrix  of  the  subject’s  velocity  may 
be  found  from  the  covariance  matrix  of  the  subject’s  state.  The  covariance  matrix  of  the 
subject’s  state  is  given  by  Equation  (2.52).  This  covariance  matrix  is  expanded  yielding 

var(p.v)  cov(pv, px)  ...  cov(p,,p*) 
cov(p*,py)  var(pv)  ...  cov(p,,py) 

— 

cov(px,  p,)  .  var(p.r) 

The  upper  left  3x3  matrix  contains  the  position  covariance  matrix.  The  lower  right 
3x3  matrix  contains  the  velocity  covariance  matrix.  These  covariance  matrices  may  be 
used  to  produce  an  error  ellipsoid  for  position  and  velocity  via  the  procedure  given  in 
Section  2.5.  In  the  next  section,  the  subject  is  still  modeled  as  having  an  approximately 
linear  path  over  short  amounts  of  time;  however,  the  time  delay  is  added  to  the  model. 

3.5  Time-Delay  and  Velocity  NLO 

This  section  addresses  the  construction  of  the  time  delay  and  velocity  NLO.  The 
previous  NLO  assumed  that  the  subject’s  motion  may  be  approximately  modeled  as  having 
constant  velocity  over  a  short  length  of  time.  That  model  may  be  more  realistic  for 
fast-moving  subjects.  However,  it  ignores  another  factor  in  the  relationship  between  the 
measurements  and  the  subject:  the  time  delay. 

The  sensors  modeled  in  this  research  are  passive,  and  therefore,  they  produce  LOS 
measurements  of  the  subject  from  the  emission  radiated  by  the  subject.  The  propagation 


(3.31) 
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time  between  the  subject  and  the  sensor  is  referred  to  as  the  time  delay.  The  time  delay  is 
given  the  variable  tj  throughout  this  thesis. 

The  first  step  in  developing  the  NLO  is  to  develop  a  description  of  a  sensor’s 
measurements  in  terms  of  the  subject’s  state.  Equation  (3.14)  is  the  equation  used  to 
simulate  a  time  delay  for  a  subject  state  and  sensor  position.  In  that  case,  data  is  simulated 
and  the  subject’s  path  is  known.  The  subject’s  state  is  unknown  for  the  NLO;  however, 
the  equation  still  provides  an  equation  for  tj  in  terms  of  the  subject’s  state  x.  The  goal  is 
to  describe  the  measurements  that  a  sensor  should  ideally  receive  in  terms  of  the  subject’s 
state.  This  function  is  depicted  as  £1  (x).  Since  the  time  delay  is  explicitly  described  in 
terms  of  the  subject’s  state  in  Equation  (3.14)  it  is  not  another  element  of  the  subject’s 
state. 

For  the  velocity  NLO,  the  ideal  measurement  angles  are  a  function  only  of  the 
subject’s  position  at  the  time  of  measurement  tm.  When  time  delay  is  considered,  the 
ideal  measurements  are  rather  only  a  function  of  the  subject’s  position  at  the  time  when 
the  subject  emitted  the  radiation  (tm  -  tci)  observed  by  the  sensor.  Furthermore,  it  may  be 
desirable  to  geolocate  the  subject  at  a  time  other  than  a  measurement  time. 

This  issue  is  encountered  in  the  previous  section  where  the  velocity  NLO  is  developed. 
The  time  at  which  the  subject’s  state  is  desired  is  called  the  estimation  time  te.  For  the 
velocity  NLO,  the  subject’s  state  at  the  time  of  measurement  is  described  in  terms  of  its 
state  at  the  desired  estimation  time.  The  same  approach  is  followed  here  the  only  difference 
being  that  the  subject’s  location  at  the  time  of  measurement  is  not  the  location  that  the 
sensor’s  LOS  measures. 

Figure  3.7  develops  the  solution  to  this  problem  in  the  same  manner  as  for  the  velocity 
NLO.  The  objective  is  to  develop  the  function  Q(x(te)).  Q(x(te))  describes  the  LOS 
measurements  as  a  function  of  the  subject’s  state  at  the  estimation  time.  The  sensor’s 
measurements  at  tm  are  measurements  of  the  subject’s  location  at  tm  minus  the  time  it  took 
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for  the  radiation  to  reach  the  sensor  from  the  subject  td.  Therefore,  the  measurements  at  the 
measurement  time  is  the  LOS  to  the  subject  at  tm  -  td.  Therefore,  the  measurements  angles 
6  and  (j>  may  be  described  as 


^  =  (p (tm  -  td)  -  s (t,„) 


mte)) 

arctan  2  (if/y,  i (rxj 

0(x(t}) 

arctan  2  (  ^if/2x  +  if/2,  ifrz J 

(3.32) 


where  s (tm)  is  the  sensor’s  position  at  /,„,  p  is  the  subject’s  location,  td  is  the  time  delay,  and 
if/  is  the  difference  vector  between  the  subject  and  the  sensor. 
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Figure  3.7:  Geometry  for  time  delay  and  velocity  NLO  development.  The  velocity  is  used 
to  project  the  subject  from  its  location  at  tm  to  its  location  corresponding  to  the  sensors’ 
measurements.  It  is  then  projected  to  its  position  at  the  estimation  time. 


Equation  (3.32)  gives  the  noiseless  LOS  measurement  that  should  be  produced  in 
terms  of  the  subject’s  state  at  te  -  td.  To  estimate  the  subject’s  state  at  te  via  the  NLO, 
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the  ideal  measurements  must  be  described  in  terms  of  the  subject’s  state  at  te.  Therefore,  a 
relationship  between  the  state  at  te  -  tci  and  te  must  be  found.  It  is  assumed  that  over  short 
durations,  the  subject’s  motion  is  approximately  linear,  and  so  it  has  constant  velocity,  ie 

P  (tm  td)  —  P(^e)  P  (te  (fm  td))  (3.33) 

This  equation  still  contains  an  unknown,  the  time  delay  td .  In  Section  3.2.3, 
Equation  (3.14)  is  developed  in  which  an  equation  for  the  time  delay  is  described  in  terms 
of  the  subject’s  state  at  trn.  However,  to  write  Q(-)  as  a  function  exclusively  of  x(te),  td 
must  be  described  in  terms  of  x(te).  The  solution  for  td  from  Equation  (3.14)  is  denoted  as 
function  td(x(tmj).  Given  the  linear  path  assumption  for  the  subject,  we  have 

P (fm)  =  P (te)  +  P (fm  -  te)-  (3.34) 

Therefore,  the  time  delay  function  may  be  described  in  terms  of  the  subject’s  state  at 
te  as  given  by 


td(Mtm))  =  td(y(te)  +  p (fm  -  te),  p)  (3.35) 

For  simplicity  this  result  will  be  called  td(x(te)). 

When  Equation  (3.35)  is  coupled  with  Equation  (3.32)  and  Equation  (3.33)  it  produces 
the  following  equation  for  the  ideal  measurements  exclusively  in  terms  of  the  subject’s  state 
at  te  as  given  by 


lj/  =  [p(fe)  -  P  (te  ~  tm  +  td(x(te)))\  ~  S (tm) 


0(x(tm)) 

arctan  2  (i)/y,  if/x) 

(p(x{tm)) 

arctan  2  (  f -J 

(3.36) 
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The  linearity  assumption  in  the  NLO  results  from  finding  the  Jacobian  of  this  equation. 
Because  the  time  delay  was  described  in  terms  of  the  subject’s  state  at  the  te,  the  time  delay 
is  not  another  element  of  the  state  vector.  Therefore,  the  Jacobian  is  formed  just  as  given 
by  Equation  (3.29)  for  the  velocity  NLO.  However,  the  differentiation  of  6  and  (p  in  terms 
of  the  subject’s  state  is  not  the  same  as  it  was  for  the  Velocity  NLO  because  the  time  delay 
is  included. 

Given  this  Jacobian,  the  NLO  iteratively  improves  upon  its  initial  guess  via 
Equation  (3.20).  The  only  difference  is  in  the  calculation  of  the  Jacobian.  Lurthermore,  the 
addition  of  the  time  delay  does  not  change  the  procedure  for  producing  the  error  surface 
or  subject’s  state  covariance  matrix.  One  issue  that  has  not  been  discussed  thus  far  is  the 
selection  of  measurements  which  are  used  to  estimate  the  subject’s  state. 

3.6  Measurement  Selection 

The  NLO  linearizes  an  inherently  nonlinear  problem.  The  linearization  transforms 
the  problem  into  iteratively  solving  an  overdetermined  system  of  equations  as  described  in 
Section  2.4.  Ignoring  weighting,  the  NLO  is  set  up  as 

A £ik  «  /tAx(.  (3.37) 

where  J  is  the  Jacobian  which  provides  the  linearizing  assumption  as  described  in  the 
previous  sections.  Ax^.  is  the  difference  between  the  current  best  estimate  of  the  subject’s 
state  Xf,  and  the  next  best  estimate.  Afl  is  the  difference  between  the  measured  angles  S2 
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and  the  LOS  angles  associated  with  xk.  Equation  (3.37)  is  expanded  yielding 
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where  6j  is  the  9  angle  that  was  measured  by  the  jth  sensor  and  9j(xk)  is  the  angle  that  would 
have  been  measured  in  a  noiseless  scenario  if  the  subject’s  state  really  was  xk. 

For  the  system  to  be  overdetermined,  the  dimension  of  £1  must  be  greater  than  the 
number  of  elements  in  the  subject’s  state.  Therefore,  if  there  are  n  elements  in  the 
subject’s  state,  then  2  x  m  >  n.  Therefore,  the  static  NLO  requires  that  a  minimum 
of  two  measurements  be  used.  The  velocity  NLO  and  time  delay  NLO  requires  three 
measurements  to  be  used.  However,  more  measurements  could  be  used. 

There  may  be  many  measurements  taken  over  a  long  timespan.  All  of  these 
measurements  could  potentially  be  used,  too.  For  many  types  of  subjects,  this  produces 
large  errors  because  of  the  assumptions  inherent  to  the  various  NLOs. 

The  static  NLO  assumes  that  the  subject  is  approximately  stationary  over  the  length 
of  time  that  the  measurements  used  by  the  NLO  are  produced.  This  assumption  may  be 
terribly  inaccurate  over  a  long  period  of  time,  so  the  resulting  estimate  would  be  very 
inaccurate.  If  the  subject’s  fiightpath  is  not  approximately  linear  over  the  entire  duration 
that  it  is  observed,  the  estimate  produced  by  the  velocity  NLOs  would  be  inaccurate. 
Therefore,  only  a  subset  of  the  measurements  should  be  used  to  estimate  the  subject. 
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The  process  of  selecting  the  measurements  that  should  be  used  to  geolocate  the  subject 
is  called  windowing.  This  is  performed  via  the  function  “setWindow.m”  provided  in 
Appendix  A.  There  are  several  factors  to  consider  when  deciding  which  measurements 
should  be  used  to  geolocate  the  subject  at  a  certain  time. 

Suppose  one  sensor  produces  several  LOS  measurements  in  a  short  time.  Because 
the  sensor’s  location  has  not  changed  significantly,  these  LOS  measurements  are  not  very 
different  from  one  another,  and  so  they  may  produce  an  inaccurate  geolocation  estimate. 
In  terms  of  the  confidence  surface,  these  measurements  will  be  incapable  of  narrowing  the 
error  surface  in  range.  Such  a  scenario  may  provide  a  good  estimate  in  two  dimensions,  but 
its  estimate  in  range  may  be  very  inaccurate.  This  fact  is  illustrated  in  Figure  3.8. 

Therefore,  it  is  desirable  to  use  multiple  sensors  with  different  views  of  the  subject. 
A  diversity  in  observation  directions  provides  information  on  the  subject’s  position  in  all 
dimensions.  Such  measurements  are  referred  to  as  diverse. 


Figure  3.8:  Effect  of  using  NLO  with  several  similar  LOS  measurements.  Since  there  is 
very  little  diversity  in  the  measurements  in  range,  the  range  estimate  may  be  inaccurate.  For 
this  reason,  it  is  desirable  to  use  multiple  independent  sensors  with  different  look-angles. 


While  unlikely,  a  situation  with  limited  viewing  diversity  is  possible  with  multiple 
sensors.  This  possible  scenario  is  ignored  in  the  windowing  process.  If  multiple  sensors 
are  used  to  provide  diversity,  then  multiple  measurements  from  the  same  sensor  might  also 
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be  used.  However,  measurements  that  weren’t  taken  near  the  estimation  time  may  cause 
inaccurate  geolocation  since  the  NLOs  assumptions  become  inaccurate. 

Therefore,  three  criteria  are  used  in  performing  the  windowing.  The  window  must 
provide  a  minimum  number  of  measurements,  nfi,  some  of  which  may  be  from  the  same 
sensor.  The  measurements  must  come  from  a  minimum  number  of  unique  sensors  ns,  and 
a  maximum  window  size  tw.  When  estimating  the  subject’s  state  at  a  time  te,  the  window 
selects  measurements  as  described  follows. 

The  window  looks  within  tw/2  of  the  estimation  time.  The  window  then  selects  the 
measurements  that  are  closest  to  the  estimation  time  that  are  taken  from  ns  different  sensors 
such  that  there  are  a  total  of  measurements. 

For  example,  consider  the  following  situation.  The  subject  is  being  estimated 
at  0  seconds,  and  the  requirements  are  a  minimum  of  three  measurements  and  two 
measurements  from  unique  sensors.  Suppose  sensor  one  produces  four  measurements  at  0, 
1,  2  seconds  and  then  sensor  two  produces  a  measurement  at  5  seconds.  If  the  window  looks 
within  5  seconds  of  the  estimation  time,  all  four  of  the  measurements  will  be  used  to  satisfy 
the  unique  sensors  requirement.  However,  if  the  window  requirement  was  4  seconds,  only 
the  first  sensor’s  measurements  would  be  used,  and  a  warning  message  would  be  produced. 
More  examples  are  illustrated  in  Figure  3.9  below. 

The  x-axis  represents  time.  The  red  dots  represent  that  times  at  which  the  subject’s 
state  is  estimated.  The  other  dots  represent  the  sensors  which  produced  measurements  at  the 
times  indicated  on  the  jt-axis.  The  window  for  selecting  the  measurements  to  estimate  the 
subject  at  each  point  in  time  are  shown  as  the  red-boxes.  The  parameters  for  the  windowing 
are  shown  at  the  top  of  the  figure. 

Due  to  the  timing  of  the  LOS  generation  detailed  in  Section  3.2,  the  measurement 
constraints  could  also  be  satisfied  without  violating  the  window  requirement.  This 
simplification  is  done  to  reduce  the  number  of  random  variables  affecting  the  comparisons 
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Figure  3.9:  Three  windowing  examples.  The  x-axis  is  the  temporal  axis.  The  red  squares 
on  the  lowest  line  represent  estimation  times.  Each  row  represents  measurement  times  for 
a  unique  sensor  designed  by  a  sensor  ID  number.  The  window  for  each  estimation  time  is 
shown  as  a  red  box. 


between  the  NLO  algorithms.  The  following  section  describes  how  each  algorithm’s 
performance  is  objectively  compared. 

3.7  Algorithm  Performance  Metrics 
3.7.1  Bias. 

Metrics  for  describing  the  geolocation  performance  of  the  algorithms  need  to  be 
defined.  As  discussed  in  section  2.4,  the  NLOs  estimates  should  be  unbiased  if  their  model 
is  accurate.  The  bias  is  in  the  algorithms  is  approximated  via  MCS.  The  average  estimate 
of  the  subject’s  position  or  velocity  is  compared  to  the  subject’s  true  position  or  velocity. 
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Therefore,  the  bias  ft  is  given  by 


P  =  \\E(x)  -  x||  (3.39) 

where  x  is  an  estimate  of  the  subject’s  true  state  x.  The  next  section  provides  a  metric  that 
describes  the  average  error  in  an  estimate. 

3.7.2  Mean  Absolute  Error. 

The  absolute  error  of  a  geolocation  estimate  is  given  by  ||x  -  x||  where  x  is  the  true 
subject  location,  and  x  is  an  estimate  of  x.  This  metric  provides  the  absolute  error  of  a 
measurement  only  for  a  specific  system  and  a  specific  set  of  measurements.  The  system  in 
this  context  means  the  locations  of  the  sensors,  their  measurements  and  the  subject’s  state. 
Because  the  measurements  are  random,  the  absolute  error  for  a  specific  system  is  also 
random.  Therefore,  a  useful  metric  for  the  estimate  error  of  the  algorithm  needs  to  describe 
the  absolute  error  of  the  algorithm  for  an  arbitrary  system  and  the  random  measurements. 

While  the  distribution  of  the  measurements  is  assumed  to  be  Gaussian  as  per  Chapter 
2,  the  covariance  of  that  distribution  is  unknown.  To  remove  extra  stochastic  parameters 
from  the  algorithm’s  overall  absolute  error,  the  algorithm’s  absolute  error  is  given  for  a 
fixed  measurement  covariance.  Therefore,  the  algorithm’s  absolute  error  is  described  by 
the  Mean  Absolute  Error  (MAE)  for  a  fixed  covariance  and  a  random  system  as  given  by 

£  =  E[||x-x||]  (3.40) 

where  £  is  the  algorithm’s  MAE.  The  precision  of  a  geolocation  algorithm  will  be  defined 
similarly. 

3.7.3  Precision. 

The  precision  of  a  specific  estimate  is  analytically  described  by  the  covariance  matrix 
Ex.  It  may  be  difficult  to  compare  the  precision  entailed  by  a  covariance  matrix  by  the 
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entries  of  the  matrix.  Therefore,  it  is  desirable  to  convert  the  covariance  matrix  into  a 
single  number. 

A  covariance  matrix  describes  the  shape  of  an  ellipsoid  [14].  Therefore,  the  precision 
may  be  described  the  volume  of  the  ellipsoid.  The  volume  of  an  ellipsoid  is 

4 

V  =  -nabc  (3.41) 

where  a,  b,  and  c  are  the  lengths  of  the  principal  axes  of  the  ellipsoid  [21].  The  principal 
axes  of  an  ellipsoid  are  given  by  the  eigenvectors  and  eigenvalues  of  the  covariance  matrix 
[18].  The  eigenvectors  are  the  principles  axes,  and  their  lengths  are  given  by  the  square  root 
of  the  eigenvalues  associated  with  these  eigenvectors.  (Because  covariance  matrices  are 
positive  definite,  their  eigenvalues  are  always  positive.)  The  eigenvalues  and  eigenvectors 
may  be  found  via  Singular  Value  Decomposition  (SVD)  [19]. 

The  SVD  decomposes  the  covariance  matrix  E  into  three  parts:  U,  S,  and  W.  For  a 
symmetric  matrix  such  as  a  covariance  matrix,  U  and  W  are  identical  [31].  The  columns  of 
W  are  the  eigenvectors  of  E.  S  is  a  diagonal  matrix  of  the  singular  values  of  E.  The  singular 
values  of  a  matrix  E  are  the  square  roots  of  the  eigenvalues  of  EffE  [25].  The  objective  is  to 
associate  the  singular  values  of  E  with  the  eigenvalues  of  E.  Because  covariance  matrices 
are  Hermitian  E/7E  =  E2.  Furthermore,  if  (T,v)  is  an  eigenpair  of  E  associated  with  the 
eigenvector  v  then 


Ev  =  Tv 
E2v  =  EvT 

E2v  =  T2v.  (3.42) 
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Therefore,  if  a  matrix  E  has  an  eigenpair  (A,  v),  then  E2  has  the  eigenpair  (A2,  v).  The 


singular  value  cr  of  E  associated  with  v  is  then  given  as 

cr  =  VIF  =  \A\.  (3.43) 

From  [18],  the  lengths  of  the  ellipsoid  associated  with  the  covariance  matrix  are  the  square 
roots  of  the  eigenvalues  of  the  covariance  matrix.  Therefore,  the  lengths  of  the  ellipsoid’s 
principle  axes  are  given  by  the  square  roots  of  the  singular  values  of  the  covariance  matrix. 

Because  the  volume  metric  is  variable  with  the  system  and  the  LOS  measurements,  an 
overall  precision  metric  is  defined  as  the  expected  value  of  the  confidence  surface’s  volume 
for  a  random  system  and  a  set  measurement  covariance  matrix.  Therefore,  the  precision 
metric  is  given  as 


p  =  E[V\.  (3.44) 

While  the  precision  metric  provides  a  means  of  measuring  the  precision  of  a 
geolocation  algorithm,  this  information  has  limited  utility  if  the  covariance  matrix  is  not 
accurate. 

3.7.4  Confidence  Estimation  Accuracy. 

A  geolocation  algorithm  may  overestimate  or  underestimate  the  confidence  in  its 
estimate  of  the  subject’s  state.  Therefore,  it  is  important  to  know  whether  or  not  an 
algorithm  accurately  represents  its  confidence  estimate.  The  accuracy  of  a  geolocation 
algorithm’s  confidence  estimates  may  be  calculated  via  the  unitless  Normalized  (state) 
Estimation  Error  Squared  NEES.  The  NEES  is  given  in  [1].  The  NEES  for  a  given 
measurement  is  defined  as 


e  =  (X  -  X)TEx(X  -  X). 


(3.45) 
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If  the  geolocation  algorithm  accurately  estimates  the  confidence  in  the  subject’s  state 
(the  covariance  matrix),  then  E  [e]  will  equal  the  degrees  of  freedom  in  the  subject’s  state. 
A  larger  value  indicates  that  the  algorithm  produces  overly  confident  estimates.  Smaller 
values  mean  that  the  algorithm  underestimates  its  confidence  in  the  subject’s  state. 

This  chapter  described  the  methodology  by  which  the  research  question  is  evaluated. 
It  described  the  simulations  and  models  used  to  synthesize  data  for  the  geolocation 
algorithms.  The  next  chapter  describes  the  performances  of  the  geolocation  algorithms 
using  synthesized  data. 
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IV.  Data  &  Results 


THis  chapter  provides  the  data  and  results  for  each  of  the  algorithms  developed  in  the 
previous  chapter.  These  algorithms  are  the  static  NLO,  velocity  NLO,  and  time 
delay  and  velocity  NLO.  Each  algorithm  assumes  certain  characteristics  for  the  subject. 
For  example,  the  static  NLO  assumes  that  the  subject  is  stationary.  Furthermore,  several 
claims  were  made  regarding  NLO  performance  when  its  assumptions  are  accurate.  As 
described  in  Chapter  2,  NLO  should  produce  an  unbiased  geolocation  estimate  that  is  the 
most  likely  estimate,  and  NLO  should  accurately  represent  the  confidence  in  its  geolocation 
estimates. 

Therefore,  subsequent  sections  provide  results  with  the  objective  of  verifying  that 
the  NLOs  converge  to  this  estimate  when  their  assumptions  exactly  model  the  scenario. 
Scenarios  where  these  assumptions  are  true  are  simulated  using  the  models  described  in 
Chapter  3.  The  setup  for  these  scenarios  is  discussed  and  the  results  gathered  from  these 
scenarios  are  provided. 

After  evaluating  the  NLOs  where  their  assumptions  are  accurate,  the  NLOs  are  used  to 
geolocate  the  most  challenging  subject:  the  Space  Vehicle  (SV).  The  NLO’s  assumptions 
are  only  approximately  true  for  this  subject.  Results  on  the  relative  performance  of 
these  algorithms  on  this  subject  are  given.  The  results  of  this  section  establish  a  means 
of  determining  if  the  more  robust  NLOs  provide  improved  geolocation  accuracy,  mean 
absolute  error  and  precision.  The  first  section  of  this  chapter  discusses  the  verification 
results  for  the  static  NLO. 

4.1  Static  NLO 

The  static  NLO  assumes  that  the  subject  is  stationary.  Therefore,  simulated  scenarios 
where  the  subject  is  stationary  are  used  to  verify  the  theoretical  performance  of  the  static 
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NLO.  The  generated  scenarios  consist  of  a  random  number  of  satellites.  The  satellites’ 
positions  are  randomly  generated  such  that  their  LOS  is  not  occluded  by  the  Earth.  Each 
sensor  is  given  a  random  orbiting  radius.  The  LOS  are  then  synthesized  as  discussed  in 
Section  3.2.  A  typical  system  generated  in  this  manner  is  shown  in  Figure  4.1. 


Figure  4.1:  Typical  static  NLO  simulation.  The  magenta  points  are  satellites.  The  red 
asterisk  is  the  subject. 


Figures  4.2  and  4.3  provide  a  visual  comparison  of  the  error  surfaces  provided  by  the 
NLO  and  TA.  These  error  surfaces  were  generated  as  discussed  in  Sections  2.5  and  2.6  for 
one  simulated  scenario.  These  surfaces  are  shown  for  three  standard  deviations  of  error. 
Thus  there  was  a  97%  probability  that  the  subject  is  located  inside  of  the  error  surfaces. 

The  top  portion  of  Figure  4.2  shows  the  size  of  the  error  surface  for  the  triangulation 
method.  The  Monte  Carlo  error  surfaces  were  generated  using  a  MCS  with  10,000  sample 
points.  Compare  the  TA’s  error  surface  in  the  bottom  plot  to  the  top  plot.  The  top  plot 
shows  the  error  surface  for  the  NLO.  Observe  that  the  NLO’s  error  surface  is  much  smaller 
than  the  TA’s  error  surface.  The  smaller  size  of  the  NLO’s  error  surface  means  that  the 
NLO  provides  tighter  bounds  (better  confidence)  in  its  estimate  of  the  subject’s  position. 

Furthermore,  the  top  portion  of  Figure  4.2  illustrates  another  claim  about  the  NLO: 
that  the  calculated  error  surface  matches  its  true  error.  In  Section  2.5,  a  method  was 
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Figure  4.2:  Error  Surfaces  for  the  Static  NLO  and  the  TA.  The  NLO’s  surface  (top)  is 
smaller  (so  the  NLO  provides  better  confidence).  The  NLO’s  calculated  error  surface 
(shown  in  red  matches)  its  error  surface  produced  using  MCS  (shown  in  teal). 
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provided  for  producing  the  NLO’s  error  surface  without  using  MCS.  This  method  allowed 
the  covariance  matrix  of  the  subject’s  state  to  be  calculated  rather  than  estimated  by  Monte 
Carlo  Simulations.  In  this  figure,  the  red  surface  is  the  NLO’s  calculated  error  surface  and 
the  teal  error  surface  is  produced  by  a  MCS.  It  is  visually  clear  that  the  calculated  error 
surface  and  Monte  Carlo  error  surface  are  nearly  the  same.  This  similarity  indicates  that 
the  NLO  accurately  estimates  its  confidence.  A  comparison  of  the  relative  sizes  of  the  error 
surfaces  is  shown  in  Figure  4.3. 

In  Figure  4.3,  the  TA  and  NLO’s  error  surfaces  are  centered  on  their  estimated  subject 
position.  The  true  subject  position  is  shown  in  white.  The  NLO’s  error  surface  is  red.  The 
TA’s  error  surface  is  blue.  It  is  clear  from  this  figure  that  the  NLO’s  error  surface  is  much 
smaller  than  that  of  the  TA.  From  Section  3.7,  the  size  of  the  error  surfaces  corresponds 
to  the  precision  of  the  algorithms.  Therefore,  the  NLO  provides  better  precision  (or 
confidence)  in  this  scenario.  The  NLO  is  also  centered  much  closer  to  the  true  subject 
position  than  the  TA.  Thus,  the  NLO  is  also  more  accurate  in  this  scenario. 

This  figure  provided  a  visual  illustration  of  the  performance  error  surface  for  the  NLO 
and  TA  in  a  typical  static  scenario.  Seven  thousand  of  these  static  scenarios  were  simulated 
to  produce  statistics  on  the  performance  of  the  static  NLO  and  the  TA.  Between  3  and 
7  satellites  were  simulated  in  the  scenarios.  The  performance  of  the  TA  and  NLO  will  be 
estimated  from  these  simulations.  In  Section  3.7,  four  statistics  were  discussed:  bias.  Mean 
Absolute  Error  (MAE),  precision,  and  the  confidence  estimation  accuracy. 

The  MAE  is  the  average  distance  between  the  estimates  of  the  subject’s  position  and  its 
true  position.  The  precision  is  described  as  the  average  volume  of  the  ellipsoid  associated 
with  some  standard  deviation  of  error.  These  values  are  given  as  a  function  of  the  number 
of  sensors  is  given  in  Figure  4.4.  It  may  not  be  easy  to  relate  the  volume  of  the  ellipse  to 
its  size.  Therefore,  to  provide  better  intuition  into  the  ellipsoid  sizes  associated  with  these 
volumes,  this  figure  also  provides  the  length  of  the  semi-major  axis. 
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Figure  4.3:  Combined  static  NLO  and  TA  error  surfaces.  The  true  subject  location  is  shown 
in  black.  The  static  NLO’s  error  surface  is  shown  in  red.  The  TA’s  error  surface  is  largest 
and  is  shown  in  blue. 


The  accuracy  of  the  error  surface  estimates  still  needs  to  be  discussed.  The  estimated 
error  surface  accuracy  is  described  by  the  NEES.  This  metric  does  not  depend  on  the 
number  of  sensors  in  the  system,  so  it  is  not  shown  as  a  function  of  the  sensor  count.  The 
means  of  the  NEES  for  each  of  these  algorithms  are  given  by 

TA  MCS  NLO  MCS  NLO 
fie  3.04  3.04  3.02 

As  discussed  in  Section  3.7,  if  the  mean  of  the  NEES  ji(  converges  to  the  degrees  of  freedom 
in  the  subject  state,  then  the  algorithm  accurately  estimates  its  error  in  the  subject’s  state. 
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Figure  4.4:  Static  NLO  and  TA  MAE  and  precision.  The  TA,  NLO,  and  MCS  of  the  NLO 
are  shown  in  these  plots.  The  MCS  of  the  NLO  tracks  the  NLO  so  closely  that  it  cannot  be 
seen  in  these  plots. 
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In  this  case,  the  subject  has  three  degrees  of  freedom:  its  x,  y,  and  z  coordinates.  Thus, 
the  MCS  accurately  estimates  the  TA’s  error  surface.  The  NLO  MCS  accurately  estimates 
its  error  surface,  and  the  error  surface  calculated  by  the  NLO  matches  its  MCS.  Thus,  the 
NLO  accurately  calculates  its  error  surface  in  these  static  scenarios. 

In  Section  2.4,  it  was  claimed  that  the  NLO  produces  unbiased  estimates  of  the 
subject’s  state.  Thus,  the  NLO’s  average  estimate  of  the  subject’s  state  should  equal  the 
subject’s  true  state.  This  is  addressed  in  Figure  4.5. 


Figure  4.5:  NLO  is  an  unbiased  estimator.  LOS  measurements  about  the  true  LOS  are 
generated  via  MCS.  The  average  of  the  resulting  estimates  of  the  subject’s  location  is 
compared  to  the  true  location.  The  MCS  is  performed  for  different  numbers  of  samples. 


This  figure  shows  the  convergence  of  MCS  for  the  NLO.  Noisy  LOS  measurements 
centered  on  the  true  LOS  were  simulated.  The  NLO  was  used  with  these  measurements. 
The  average  of  the  NLO’s  estimate  of  the  subject’s  position  was  calculated  as  a  function  of 
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the  number  of  Monte  Carlo  Samples.  This  figure  illustrates  that  as  more  samples  are  used, 
the  average  NLO  estimate  approaches  the  subject’s  true  position. 

This  section  provides  and  discusses  the  results  for  the  static  NLO  algorithm.  The 
algorithm  was  tested  in  many  scenarios  where  all  of  the  Static  NLO’s  assumptions  are  true: 
where  the  subject  is  stationary.  Under  these  conditions  the  simulated  performance  should 
match  the  theoretical  performance.  A  similar  discussion  is  provided  for  the  velocity  NLO 
in  the  next  section. 

4.2  Velocity  NLO 

As  in  the  previous  section,  the  objective  of  this  section  is  to  demonstrate  that  the 
simulated  performance  of  the  velocity  NLO  matches  the  theoretical  NLO.  The  velocity 
NLO  is  based  on  the  assumption  that  the  subject’s  motion  is  approximately  linear  over  short 
durations  of  time.  Furthermore,  sensor  measurements  are  not  necessarily  simultaneous. 
These  scenarios  are  simulated  as  described  in  Section  3.4. 

A  subject  is  simulated  with  a  linear  flightpath  and  a  constant  velocity.  Realistic  sensors 
and  their  noisy  measurements  are  also  simulated.  The  velocity  NLO  is  used  with  this  data, 
and  statistics  are  produced. 

As  described  in  Section  3.6,  when  sensor  measurements  are  not  simultaneous,  a  time  at 
which  to  estimate  the  subject’s  state  must  be  chosen.  Furthermore,  selecting  measurements 
from  which  to  estimate  the  subject’s  state  are  selecting  via  the  measurement  window.  The 
measurement  window  has  three  settings.  The  settings  used  to  produce  the  data  for  the 
Velocity  NLO  are  given  below. 

A  total  of  five  satellites  were  generated.  A  minimum  of  seven  measurements  were  used 
to  perform  the  NLO.  These  seven  measurements  must  come  from  at  least  three  independent 
satellites,  and  the  time  window  was  set  at  20  seconds.  In  addition  to  using  the  velocity 
NLO,  the  static  NLO  and  TA  were  also  used  to  geolocate  the  constant  velocity  subject. 
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The  static  NLO  and  TA  do  not  account  for  velocity,  so  the  MAE  of  their  geolocation 
estimates  should  degrade  as  the  subject’s  velocity  increases.  Because  only  velocity  has 
been  added  to  the  system  model,  these  statistics  illustrate  the  added  MAE  due  to  velocity 
in  the  static  NLO  and  TA. 

To  observe  these  effects  as  a  function  of  the  subject  velocity,  the  subject’s  velocity  was 
sampled.  A  constant  velocity  subject  was  simulated  with  magnitudes  of  10,  50,  100,  500, 
1000,  3000,  and  7000  meters  per  second.  At  each  velocity,  twenty  scenarios  were  produced 
from  which  more  than  10,000  state  estimates  were  recorded  to  generate  statistics. 

Figure  4.6  shows  the  MAE  in  position  estimates  for  the  three  geolocation  algorithms 
versus  subject  velocity.  Note  that  the  constant  velocity  NLO’s  estimates  maintain  the  same 
MAE  irregardless  of  increased  subject  velocity.  This  is  expected  since  the  NLO  optimizes 
over  velocity.  Further  observe  that  the  static  NLO  provides  the  smallest  MAE  for  very 
slow  subjects;  however,  its  MAE  degrades  the  fastest.  Once  the  subject’s  velocity  is  greater 
than  1000  meters  per  second,  the  static  NLO  appears  to  diverge.  This  means  that  the  static 
NLO’s  assumption  that  the  subject  is  approximately  stationary  is  so  inaccurate  that  it  no 
longer  converges.  The  TA’s  estimates  are  more  accurate  than  the  static  NLO  at  faster 
velocities,  but  it  performs  worse  than  the  velocity  NLO. 

The  MAE  of  the  velocity  estimates  is  similar  to  the  error  for  the  position  estimates. 
This  is  shown  in  Figure  4.7.  As  expected,  the  velocity  NLO  provides  constant  MAE  in 
its  velocity  estimates.  Observe  the  same  divergent  behavior  in  static  NLO,  and  the  TA’s 
decreased  MAE.  Note  that  neither  the  static  NLO  nor  the  TA  innately  calculate  velocity. 
Their  velocity  estimates  were  calculated  by  interpolation.  A  remaining  question  is  the  mean 
accuracy  of  the  error  surface  estimates  for  the  NLOs. 

Figure  4.8  shows  the  NEES  for  the  static  NLO  as  a  function  of  subject  velocity.  The 
NEES  starts  at  3.1  (slightly  over  confident)  and  ends  with  a  value  near  108.  For  such  large 
velocities,  the  static  NLO  is  drastically  over  confident.  The  NEES  of  the  velocity  NLO  was 
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Position  Estimation  Error 


Figure  4.6:  Mean  absolute  error  in  the  estimated  position  for  a  constant  velocity  subject 
without  time  delay.  The  velocity  NLO  has  a  constant  mean  absolute  error  as  expected.  The 
TA  and  the  static  NLO’s  mean  absolute  error  degrades  with  increased  velocity. 


estimated  to  be  six.  Note  that  since  the  velocity  NLO  optimizes  over  position  and  velocity, 
it  estimates  six  parameters  of  the  subject’s  state.  Therefore,  an  NEES  of  six  is  ideal.  The 
velocity  NLO  accurately  calculates  its  error  surface  for  all  subject  velocities. 

The  velocity  NLO  should  also  be  an  unbiased  estimator  meaning  that  its  geolocation 
estimates  should  be  unbiased.  Furthermore,  because  the  static  NLO  does  not  include 
subject  motion  in  its  model,  its  geolocation  estimates  may  be  biased.  Therefore,  MCS 
were  used  where  the  LOS  measurements  of  the  subject  were  simulated  about  the  true  LOS. 
The  static  and  velocity  NLOs  were  used  with  this  data  to  geolocate  the  subject.  The  result 
of  this  MCS  is  shown  in  Figure  4.9. 
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Velocity  Estimation  Error 


Figure  4.7:  Mean  absolute  error  in  velocity  estimates  for  the  velocity  NLO,  static  NLO,  and 
TA.  The  velocity  NLO  performs  as  expected.  The  error  in  the  other  algorithms’  geolocation 
estimates  increases  with  the  subject’s  velocity. 


The  results  provided  in  this  section  describe  the  performance  of  the  TA  and  static 
and  velocity  NLOs.  The  velocity  NLO’s  simulated  performance  matches  its  theoretical 
performance.  The  next  section  describes  a  similar  comparison  among  the  algorithms  when 
the  time  delay  is  included. 

4.3  Time-Delay  and  Velocity  NLO 

This  section  examines  the  performance  of  the  time  delay  and  velocity  NLO.  For 
brevity  in  this  section,  the  time  delay  and  velocity  NLO  will  be  referred  to  as  NLOtd- 
As  in  the  previous  sections,  the  primary  objective  of  this  section  is  to  determine  if  the 
simulated  performance  of  the  NLOtd  matches  its  theoretical  performance.  Furthermore, 
the  previous  algorithms  do  not  account  for  time  delay,  so  this  section  illustrates  the  added 
error  in  the  previous  algorithms  due  to  the  time  delay. 
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Figure  4.8:  The  static  NLO’s  NEES  for  the  constant  velocity  subject  without  time  delay. 
The  NEES  grows  as  the  subject’s  speed  increases.  This  means  that  the  static  NLO 
drastically  over  estimates  its  confidence  in  its  estimates. 


The  simulated  subject  in  this  section  is  identical  to  the  linear  subject  in  the  previous 
section,  except  that  the  time  delay  affects  the  LOS  measurements.  These  scenarios  consist 
of  a  single  subject  moving  with  constant  velocity.  It  is  observed  by  five  sensors.  To  estimate 
the  subject,  a  minimum  of  seven  measurements  are  used  from  at  least  three  unique  sensors. 
All  of  these  measurements  must  be  found  within  20  seconds  of  the  estimation  time. 

If  unaccounted  for,  the  time  delay  creates  a  bias  in  the  LOS  measurements.  This  bias 
increases  with  the  subject’s  velocity.  Bias  is  induced  into  the  measurements  because  the 
subject  moves  some  distance  during  the  time  that  it  takes  for  its  emission  to  be  observed  by 
the  sensors.  The  faster  the  subject  moves,  the  further  it  has  moved  over  this  time.  Therefore, 
a  LOS  measurement  produced  at  a  certain  time  will  not  correspond  to  a  measurement  of  the 
subject  at  that  time.  Thus  the  LOS  measurements  lag  the  position  of  the  subject  creating  the 
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NLO  MCS  Convergence 


Figure  4.9:  The  velocity  NLO  is  an  unbiased  estimator  for  a  linear  subject  without  time 
delay. 


bias.  The  following  result  elucidates  the  magnitude  of  this  error  in  estimating  the  subject’s 
location  as  a  function  of  velocity. 

Subjects  were  simulated  with  velocities  of  10,  50,  100,  500,  1,000,  3,000,  and  7,000 
meters  per  second.  Statistics  at  each  of  these  velocities  are  based  on  more  than  10,000 
geolocation  samples.  The  TA,  static  and  velocity  NLOs,  and  NLOtd  were  used  to  estimate 
the  subject. 

Figure  4.10  shows  the  MAE  of  the  position  estimates  for  each  geolocation  algorithm 
as  a  function  of  the  subject’s  velocity.  TheMAE  of  the  NLOtd  should  be  unaffected  by 
velocity  because  it  accounts  for  the  time  delay.  The  other  three  algorithms  do  not,  and  so 
their  performance  should  degrade  as  the  subject’s  velocity  increases.  Once  again,  the  static 
NLO  performs  the  best  for  slower  subjects;  however,  its  performance  degrades  quickly  as 
the  subject’s  speed  increases.  The  velocity  NLO  performs  nearly  as  well  as  the  NLOTd 
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for  all  but  the  fastest  subjects.  For  the  fastest  subject,  the  NLOTD  provides  an  MAE  of  147 
meters  versus  654  meters  for  the  velocity  NLO  and  4,691  meters  for  the  TA.  Similar  results 
are  shown  for  the  velocity  estimates. 


Figure  4.10:  The  MAE  of  each  geolocation  algorithm’s  position  estimates.  The  time 
delay  and  velocity  NLO  provides  constant  performance.  The  other  algorithms’  accuracies 
degrade  with  increased  subject  velocity. 


These  effects  are  also  seen  in  the  velocity  plot  shown  in  Figure  4.11.  As  before, 
the  performance  of  the  NLOtd  is  unaffected  by  the  speed  of  the  subject.  All  of  the 
other  algorithms  produce  less  accurate  estimates  as  the  subject’s  speed  increases.  The 
velocity  NLO  matches  the  NLOTd  MAE  much  better  than  it  matches  position.  Now  that 
the  geolocation  MAE  has  been  addressed  for  all  three  NLOs,  their  NEES  is  given. 

Figure  4.12  shows  the  NEES  of  the  static  and  velocity  NLO  as  functions  of  the 
subject’s  velocity.  Note  that  the  static  NLO  and  velocity  NLOs  ideally  have  NEES  values 
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Figure  4.11:  The  MAE  of  each  geolocation  algorithm’s  velocity  estimates.  The  NLOs 
that  include  velocity  and/or  the  time  delay  provide  a  nearly  constant  MAE.  The  other 
algorithms’  MAEs  degrade  with  increased  subject  speed. 


of  3  and  6,  respectively.  The  static  NLO’s  NEES  begins  to  grow  for  subjects  with  speeds 
of  10  m/s.  The  velocity  NLO  provides  accurate  confidence  over  a  larger  range  of  subject 
velocities  than  the  static  NLO;  however,  its  NEES  also  grows  with  the  subject’s  velocity. 
This  indicates  that  the  confidence  estimates  given  by  the  static  and  velocity  NLOs  are 
overconfident.  As  expected,  the  NLOTD  was  estimated  to  have  a  NEES  of  6  for  all  subject 
velocities. 

Because  the  Velocity  and  Time  Delay  NLO  includes  time  delay  and  the  subject’s 
velocity  in  its  model,  its  estimates  should  be  unbiased.  Furthermore,  the  other  NLOs  do 
not  include  the  time  delay  in  their  models,  so  their  geolocation  estimates  may  be  biased. 
Therefore,  a  MCS  was  performed  where  LOS  measurements  of  the  subject  were  simulated 
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Figure  4.12:  The  NEES  of  the  static  and  velocity  NLOs  for  a  linear  subject  with  time  delay 
included. 


about  the  true  LOS.  The  NLOs  were  used  with  this  data  to  geolocate  the  subject.  The  result 
of  this  MCS  is  shown  in  Figure  4.13. 

This  section  addressed  the  performance  of  the  NLOtd  and  the  other  geolocation 
algorithms  on  a  subject  that  moves  with  constant  velocity.  This  subject  was  designed  to 
match  the  assumptions  inherent  to  the  NLOTD.  The  next  section  describes  the  performance 
of  the  geolocation  algorithms  on  a  most  challenging  subject. 

4.4  Relative  NLO  Performance 

In  Section  3.1.2,  an  SV  is  described  as  the  most  challenging  subject  for  geolocation. 
This  subject  is  anticipated  to  maximally  induces  the  errors  in  geolocation  due  to  the 
subject’s  path  and  time  delay.  The  previous  three  sections  describe  the  performance  of  the 
NLOs  on  simplistic  subjects  designed  to  verify  the  anticipated  performance  of  the  NLOs. 
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NLO  MCS  Convergence 


— * —  NLOtd  Position 
— -* —  NLOTD  Velocity 

— ♦ — Velocity  NLO  Position 
— -  Velocity  NLO  Velocity 
— • —  Static  NLO 


MCS  Samples 


Figure  4.13:  The  velocity  and  time  delay  NLO  is  unbiased.  The  other  geolocation 
algorithms  produced  biased  estimates. 


This  section  provides  results  that  will  be  used  to  compare  the  NLOs  and  TA.  A  typical 
simulated  SV  scenario  is  shown  in  Figure  4.14. 

Because  this  section  gives  results  for  the  most  challenging  subject,  these  results  give 
the  worst  geolocation  performance  that  are  are  expected  for  any  subject.  As  discussed  in 
Section  3.2,  the  LOS  measurements  were  arbitrarily  given  a  standard  deviation  of  error 
of  5  micro  radians.  The  performance  is  dependent  on  the  accuracy  of  the  sensors.  For 
example,  the  accuracy  and  precision  of  the  geolocation  algorithms  will  decrease  if  the 
standard  deviation  of  error  in  the  measurements  increases. 

To  provide  the  statistics  shown  in  this  section,  an  SV’s  flightpath  is  simulated  as 
described  in  Section  3.1.2.  Results  are  shown  for  3,  4,  5,  6,  7,  and  8  sensors.  A  minimum  of 
seven  measurements  were  used  from  at  least  3  independent  sensors  within  a  20  seconds  of 
the  estimation  time.  Twenty  different  SV  scenarios  were  simulated  for  each  sensor  count. 
Each  scenario  typically  provides  more  than  1,500  sensors  measurements.  The  subject  is 
estimated  at  each  measurement  time.  Thus,  these  statistics  are  based  on  approximately 
30,000  geolocation  estimates. 
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Figure  4.14:  A  typical  SV  simulation.  The  blue  points  represent  satellites.  Their  measured 
LOS  are  shown  as  the  magenta  lines.  The  yellow  line  is  the  subject’s  path,  and  the  point  in 
red  is  the  subject’s  current  position. 


Figure  4.15  shows  the  MAE  of  the  geolocation  algorithms  as  a  function  of  the  number 
of  sensors.  The  staticNLO  seems  to  diverge  regardless  of  the  number  of  sensors  present. 
The  TA’s  estimates  have  a  smaller  MAE  than  the  static  NLOs,  and  its  MAE  decreases 
with  an  increased  number  of  sensors.  The  velocity  NLOs  estimates  are  nearly  ten  times  as 
accurate  as  the  TAs.  Its  MAE  also  slowly  increases  with  more  sensors.  The  most  advanced 
NLO,  the  time  delay  and  velocity  NLO  produces  the  most  accurate  estimates.  The  MAE 
of  its  estimates  also  increase  when  using  more  sensors. 

The  same  effects  are  seen  in  the  velocity  estimates  provided  by  the  NLOs.  The  velocity 
estimates  are  shown  in  Figure  4.16.  As  before,  the  static  NLO  appears  to  diverge  except 
when  there  are  more  than  six  sensors  present.  The  TA’s  performance  is  constant  with  an 
increased  number  of  sensors.  The  velocity  NLO  and  time  delay  and  velocity  NLO  provide 
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Figure  4.15:  The  MAE  of  each  geolocation  algorithms’  position  estimates. 


similar  performance.  However,  the  time  delay  and  velocity  NLO  provides  the  best  estimate 
of  the  subject’s  velocity  for  each  number  of  sensors. 

Figure  4.17  shows  the  mean  of  the  NEES  for  each  geolocation  algorithm.  The 
static  NLO’s  NEES  is  much  greater  than  its  ideal  value  of  three.  This  indicates  that  its 
confidence  estimates  are  overconfident.  The  velocity  NLO’s  confidence  estimates  are 
also  overconfident  irrespective  of  the  number  of  sensors  present.  The  time  delay  NLO’s 
confidence  estimates  are  also  overconfident.  However,  its  error  surface  estimates  become 
much  more  accurate  as  the  number  of  sensors  increases.  With  eight  sensors,  the  time  delay 
and  velocity  NLO  has  an  NEES  of  6.1.  Ideally  it  should  have  a  value  of  6,  so  it  is  slightly 
overconfident. 

These  results  were  all  produced  using  an  arbitrarily  chosen  sensor  confidence.  Each 
sensor’s  confidence  is  described  by  the  standard  deviation  of  error  in  its  distributions  for 
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Figure  4.16:  The  MAE  of  each  geolocation  algorithms’  velocity  estimates. 


6  and  (p.  The  previously  results  used  a  value  of  5/uR.  Figure  4.18  gives  the  performance 
of  the  time  delay  and  velocity  NLO  as  a  function  of  the  standard  deviation  of  error  in  the 
sensors.  Eight  sensors  were  used  with  the  same  window  settings  as  the  previous  results. 
This  data  describes  how  the  results  scale  with  the  sensor  confidence.  These  results  are  only 
shown  for  the  time  delay  and  velocity  NLO  since  it  provides  the  most  accurate  geolocation 
estimates  for  the  SV. 

The  previous  discussion  concludes  the  results  for  this  research.  This  chapter  first 
verified  that  each  NLO  matched  its  theoretical  performance  in  simulations.  It  also  gave 
information  on  how  velocity  and  time  delay  increase  the  error  in  geolocation  when  these 
are  not  included  in  the  NLO  model.  Lastly,  results  were  given  detailing  the  performance  of 
each  geolocation  algorithm  on  the  most  challenging  subject.  The  final  chapter  will  provide 
conclusions  reached  from  this  data,  and  offer  future  research  recommendations. 


91 


Mean  of  NEES 


Figure  4.17:  The  mean  of  the  NEES  for  each  geolocation  algorithm’s  estimates. 
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Figure  4.18:  The  time  delay  and  velocity  NLO’s  performance  as  a  function  of  sensor 
confidence. 
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V.  Conclusions  &  Discussion 


In  this  final  chapter,  conclusions  derived  from  the  results  given  in  the  last  chapter  are 
discussed.  The  first  section  of  this  chapter  describes  conclusions  regarding  theoretical 
claims  about  the  geolocation  algorithms.  The  next  section  describes  conclusions  regarding 
the  magnitude  of  the  errors  resulting  from  geometry,  velocity,  and  time  delay.  It  also 
discusses  the  performances  of  the  algorithms  when  applied  to  the  most  challenging  subject. 
Lastly,  recommendations  and  discussion  are  given  for  future  work. 

5.1  Theoretical  Verification 

In  Chapter  2  several  claims  were  made  regarding  the  geolocation  algorithms.  These 
theoretical  claims  were  tested  via  simulations.  The  motivation  for  investigating  Non- 
Linear  Otpimization  (NLO)  was  to  resolve  anticipated  shortcomings  in  the  Triangulation 
Algorithm  (TA).  The  TA  does  not  account  for  system  geometry,  subject  motion,  or  the  time 
delay. 

In  Section  4.1,  a  scenario  was  simulated  with  a  stationary  subject.  The  TA  and  static 
NLO  (which  includes  system  geometry  in  its  model)  were  used  to  estimate  the  subject’s 
location.  Figure  4.4  gave  the  accuracy  of  the  geolocation  and  confidence  estimates  for  the 
TA  and  static  NLO.  It  was  shown  in  this  section,  that  the  NLO’s  geolocation  estimates 
have  half  the  MAE  provided  by  the  TA’s  estimates.  Furthermore,  the  NLO  provided  more 
confidence  in  its  estimates  (its  error  surfaces  were  smaller).  Both  algorithms’  metrics 
improved  with  more  sensors.  These  results  demonstrate  that  the  NLO  outperforms  the 
TA  by  accounting  for  system  geometry. 

Furthermore,  the  static  NLO  should  be  an  unbiased  estimator  for  a  static  subject.  Thus, 
the  static  NLO’s  average  estimate  is  the  true  subject  location.  This  claim  was  tested  in 
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Figure  4.5.  The  figure  demonstrated  that  the  static  NLO  is  an  unbiased  estimator  in  the 
static  case. 

The  static  NLO  also  accurately  calculates  its  error  surface  (the  confidence  in  its 
estimates)  as  illustrated  in  Figure  4.2.  This  was  also  demonstrated  by  the  mean  of  its 
NEES.  This  figure  demonstrated  that  the  static  NLO’s  calculated  error  surface  matches  the 
error  surface  produced  by  MCS.  However,  these  results  should  occur  only  when  the  static 
NLO’s  assumptions  are  exactly  true.  The  static  NLO  should  not  produce  ideal  results  when 
other  factors  are  included  such  as  the  subject’s  motion  and  time  delay. 

The  velocity  NLO  was  tested  in  Section  4.2.  In  this  section,  a  subject  was  simulated 
with  constant  velocity.  The  static  NLO  and  TA  assume  that  the  subject  is  stationary  while 
the  LOS  measurements  needed  to  estimate  it  are  produced.  As  this  approximation  becomes 
less  accurate,  the  static  NLO  and  TAs  geolocation  estimates  have  an  increased  MAE. 

This  hypothesis  was  tested  in  Ligures  4.6,  4.7,  and  4.8  by  implementing  the  algorithms 
on  subjects  moving  at  different  speeds.  The  results  showed  that  the  static  NLO  out¬ 
performed  the  TA  and  velocity  NLO  for  subjects  with  speeds  less  than  100  m/s.  The  static 
NLO  out-performs  the  velocity  NLO  for  slow  subjects  because  the  velocity  NLO  matches 
its  velocity  estimate  to  the  random  noise. 

As  the  subject’s  velocity  increases,  the  static  NLO’s  estimates  have  the  largest  MAE. 
This  degradation  occurs  because  the  static  NLO’s  assumptions  become  less  accurate  as  the 
subject’s  speed  increases.  It  appears  that  the  static  NLO  may  have  failed  to  converge  for 
subjects  with  speeds  greater  than  1,000  m/s.  As  anticipated,  the  TA’s  MAE  also  decreased 
with  speed,  although  at  a  slower  rate  than  the  static  NLO. 

The  velocity  NLO  accounts  for  velocity,  so  the  subject’s  velocity  is  innately  estimated. 
Lor  subject’s  with  speeds  greater  than  100  m/s,  the  velocity  NLO  provides  the  smallest 
MAE  in  its  geolocation  estimates.  Its  MAE  is  constant  for  all  subject  velocities.  The 
velocity  NLO  should  also  be  an  unbiased  estimator.  Ligure  4.9  demonstrated  that  the 
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velocity  NLO  is  an  unbiased  estimator.  Furthermore,  the  figure  showed  that  the  static 
NLO  produces  biased  geolocation  estimates  for  a  moving  subject. 

The  velocity  NLO  should  also  accurately  estimates  the  confidence  in  its  estimates.  The 
mean  of  its  NEES  was  calculated  to  be  six.  Because  the  velocity  NLO  estimates  the  position 
and  velocity  of  the  subject  in  the  x,  y,  and  z  directions,  it  estimates  six  parameters  of  the 
subject’s  state.  Therefore,  the  velocity  NLO  accurately  calculated  the  error  surface  (the 
confidence  in  its  geolocation  estimates)  in  this  scenario.  However,  Figure  4.8  demonstrates 
that  the  static  NLO  over  estimates  its  geolocation  confidence,  so  its  error  surfaces  are  not 
accurate.  The  static  NLOs  geolocation  estimates  are  overconfident  because  the  static  NLO 
does  not  include  velocity  in  its  model.  Similar  results  were  shown  for  the  time  delay  and 
velocity  NLO. 

Section  4.3  addressed  the  performance  of  the  time  delay  and  velocity  NLO.  A 
subject  was  simulated  with  constant  velocity  and  time  delay.  The  error  resulting  from  the 
time  delay  increases  as  the  subject’s  velocity  increases;  therefore,  a  faster  subject  should 
accentuate  this  error.  The  accuracy  results  for  each  algorithm  was  given  as  a  function  of 
the  subject’s  velocity  to  illustrate  this  effect.  The  TA,  static  NLO,  and  velocity  NLO  do  not 
take  the  time  delay  into  account  and  so  their  estimates  should  become  less  accurate  as  the 
subject’s  velocity  increases. 

This  effect  was  observed  in  Figures  4.10  4.11,  and  4.12.  As  anticipated,  the  time 
delay  and  velocity  NLO  provides  a  smaller  MAE  and  improved  confidence  estimates  as  the 
subject’s  velocity  increases  compared  to  the  other  algorithms.  The  static  NLO  performs 
the  best  for  subject’s  with  speeds  less  than  100  m/s.  However,  the  static  NLO’s  confidence 
estimates  are  overconfident.  The  velocity  NLO  also  overestimates  its  confidence  by  a  lesser 
amount.  For  subjects  with  speeds  less  than  100  m/s,  the  velocity  NLO  accurately  estimates 
its  error  surfaces  since  its  NEES  is  approximately  six  over  this  region.  The  time  delay  and 
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velocity  NLO  accurately  estimates  its  geolocation  confidence  irrespective  of  the  subject’s 
speed.  The  accuracy  of  its  position  and  velocity  estimates  is  also  constant. 

Furthermore,  the  time  delay  and  velocity  NLO  should  be  an  unbiased  estimator  in 
this  case.  This  was  addressed  by  Figure  4.12.  The  figure  demonstrates  the  time  delay  and 
velocity  NLO  is  an  unbiased  estimator.  The  other  NLOs  were  biased  because  they  do  not 
include  time  delay  in  their  models. 

The  aforementioned  results  demonstrated  the  error  resulting  from  the  system 
geometry,  subject  motion,  and  time  delay.  For  example,  the  faster  the  subject,  the  more 
error  is  produced  in  geolocation  estimates.  This  may  give  validity  to  the  assumption  that 
an  SV  is  the  most  challenging  subject  for  geolocation.  An  SV  is  the  fastest  subject  that 
might  be  geolocated  and  its  motion  is  complex  (meaning  it  accelerates  and  has  other  higher 
order  moments).  Therefore,  it  could  be  reasonably  expected  to  produce  the  most  error  in 
geolocation  estimates.  The  application  of  these  geolocation  algorithms  to  such  a  subject  is 
discussed  in  the  next  section. 

5.2  Algorithm  Comparison 

The  ultimate  test  of  the  geolocation  algorithms  considered  here  was  to  apply  them 
against  the  most  challenging  subject:  a  Space  Vehicle  (SV).  Section  4.4  describes  the 
results  of  applying  these  algorithms  to  geolocate  simulated  SVs.  Figures  4.15  and  4.16 
demonstrated  that  the  static  NLO  does  not  converge  to  the  estimate  nearest  to  subject’s 
position  for  SVs.  This  divergence  indicates  that  its  assumptions  are  too  inaccurate  to 
geolocate  such  a  subject.  The  TA’s  estimates  have  a  larger  MAE  than  the  velocity  and  time 
delay  and  velocity  NLOs.  The  TA’s  MAE  slightly  improves  with  an  increased  number  of 
sensors. 

The  MAE  of  the  velocity  NLO’s  estimates  is  nearly  one  hundred  times  smaller  than 
the  MAE  for  the  TA’s  estimates.  The  time  delay  and  velocity  NLO  provides  the  smallest 
MAE  in  position  and  velocity.  These  results  were  produced  for  a  sensor  confidence  of 
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5/rR.  Figure  4.18  provides  results  for  the  time  delay  and  velocity  NLO  with  other  sensor 
confidences.  As  anticipated,  more  accurate  sensors  result  in  smaller  MAE  in  geolocation 
estimates  and  improved  confidence.  Thus  changing  sensor  confidence  results  in  a  vertical 
shift  in  the  geolocation  accuracy  plots  given  in  this  research. 

Perhaps  most  important,  the  time  delay  and  velocity  NLO’s  NEES  nearly  converges 
to  its  ideal  value  for  a  sufficient  number  of  sensors.  Figure  4.17  demonstrates  that  none  of 
the  other  geolocation  algorithms  converge  to  their  ideal  NEES  value.  Therefore,  the  time 
delay  and  velocity  NLO  provides  the  most  accurate  geolocation  confidence  estimates  for 
the  most  challenging  subject. 

While  the  static  NLO  could  not  provide  useful  estimates  for  the  most  challenging 
subject,  its  estimates  have  the  smallest  MAE  for  the  position  and  velocity  of  slow  subjects. 
Therefore,  if  the  type  of  subject  being  observed  is  known,  then  it  may  be  better  to  use  the 
static  NLO.  For  example,  if  the  user  knows  they  are  observing  a  static  subject,  then  the 
static  NLO  may  provide  better  results. 

The  MAE  of  the  velocity  NLO’s  estimates  is  larger  than  the  static  NLO’s  MAE  for 
low-speed  subjects  and  larger  than  the  time  delay  and  velocity  NLO’s  MAE  for  high¬ 
speed  subjects.  Therefore,  the  velocity  NLO  does  not  appear  to  be  useful.  For  unknown 
subjects,  it  may  be  best  to  use  the  time  delay  and  velocity  NLO  because  it  is  the  most 
robust  algorithm  provided  here  and  gives  the  most  accurate  confidence  estimates  of  the 
geolocation  algorithms. 

The  previous  discussion  answers  the  first  part  of  the  research  question:  which  NLO 
performs  the  best?  The  static  NLO  provides  the  best  geolocation  algorithm  for  subjects 
with  speed  less  than  10  m/s.  The  time  delay  and  velocity  NLO  provides  the  best  geolocation 
performance  for  general  subjects.  The  second  part  of  the  research  question  is  how 
confidence  in  estimates  might  be  visualized.  A  visualization  of  geolocation  confidence 
was  provided  by  the  error  surfaces.  This  research  also  demonstrated  that  the  error  surfaces 
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can  accurately  depict  the  confidence  in  geolocation  estimates.  However,  none  of  the 
algorithms  could  accurately  estimate  the  confidence  in  their  geolocation  estimates  for  the 
most  challenging  subject.  The  research  question  has  been  addressed;  however,  there  are 
still  remaining  problems  that  may  be  addressed  in  follow-on  research. 

5.3  Future  Work 

As  described  in  Chapter  1,  there  were  many  assumptions  and  limitations  for  this 
research.  Perhaps  the  most  significant  of  the  assumptions  is  that  the  measurements  are 
unbiased.  If  LOS  measurements  are  biased,  then  none  of  the  results  given  here  will 
hold.  In  this  case,  the  NLOs  will  all  over-estimate  their  confidences.  Depending  on  the 
magnitude  of  the  bias,  the  bias  could  be  a  larger  source  of  error  than  the  random  jitter  in 
the  measurements. 

Furthermore,  the  sensors  in  this  research  measured  the  LOS  with  a  pseudo-random 
periodicity.  Measurement  times  may  not  be  consistent  in  reality.  The  measurement  times 
were  also  assumed  to  be  perfectly  known.  In  reality,  there  may  be  uncertainty  in  the 
measurement  times.  It  may  be  valuable  to  analyze  the  effects  of  adding  more  realism  to  the 
sensor  measurement  times  on  geolocation.  Relativistic  effects  could  also  be  considered. 

The  effects  of  different  sensors  on  geolocation  have  not  been  addressed.  For  example, 
it  may  be  valuable  to  consider  the  accuracy  and  confidence  improvement  resulting  from 
adding  one  sensor  close  to  the  subject.  The  accuracy  and  confidence  in  geolocation 
estimates  could  be  calculated  as  a  function  of  sensor  distance  to  the  subject. 

Lastly,  the  effects  of  filtering  the  LOS  data  could  be  useful.  While  the  noise  in 
the  LOS  measurements  is  independent,  the  measurements  are  still  taken  of  the  same 
subject.  Therefore,  one  LOS  measurement  provides  information  on  the  next  measurement. 
A  Kalman  filter  could  be  used  to  improve  the  LOS  measurements.  Something  akin  to 
this  began  to  be  examined;  however,  it  was  mistakenly  performed.  A  discussion  of  this 
investigation  is  given  in  Appendix  C. 
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These  are  just  some  of  the  ideas  that  arose  over  the  course  of  this  research.  The 
conclusions  and  discussion  provided  here  give  closure  to  the  research  question,  and  have 
resulted  in  potential  queries  for  further  research. 
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Appendix  A:  Code 


A.l  SV  _PathGen.m 


1  function  [  sv  ]  =  sv_PathGen (  sv  ) 

2  %%  Authorship 

3  %  Stephen  Hartzell,  AFIT  -  GRA 

4  %  stephen.hartzell.ctr@afit.edu,  hartzell.stephen@gmail.com 

5  %  Created:  06-01-2011,  Modified:  10-27-2012 

6  % 

7  %  sv_PathGen  Generate  or  update  an  space  vehicle  (SV)  flightpath 

8  %  [  sv  ]  =  sv.PathGen (  sv  )  performs  one  of  two  operations 

9  %  depending  on  the  format  of  the  input  sv.  If  input  sv  is  not  a 

10  %  structure,  then  realistic  parameters  that  completely  define  . . . 

an  sv ' s 

11  %  path  are  generated.  If  input  sv  is  a  structure  with  the  required 

12  %  elements,  then  the  subject's  positions  and  velocities  at  the  ... 

times 

13  %  defined  by  sv.t  are  generated  for  an  sv's  flight  path  in  a 

14  %  2D— plane.  The  resulting  sv  is  then  rotated  according  to  a  ZXY 

15  %  rotation  and  the  new  3D  points  are  stored  in  sv.r  structure.  The 

16  %  full  sv  structure  will  have  up  to  fifteen  <lxl>  structure  . . . 

fields : 

n  %  a.  Alpha,  b.  Beta,  C_e,  Omega,  Rot_X,  Rot_Y,  Rot_Z,  t,  N,  r,  v, 

is  %  t_reentry_max ,  and  t_max 

19  % 

20  %  N  is  a  vector  pointing  from  center  of  the  Earth  towards  a  ... 

plane  which  it 

21  %  is  normal  to.  The  sv's  flight  path  will  be  generated  such  . . . 

that  the 

22  %  plane  containing  the  flight  path  is  orthogonal  to  the  plane  . . . 

which  N  is 

23  %  orthogonal  to.  N  is  useful  if  the  sv  flight  path  must  be  . . . 

located  such  that 

24  %  it  is  observable  by  a  system  of  satellites. 

25  % 

26  %  NOTE  1:  The  only  known  problem  with  this  generator  is  that  . . . 

for  large 

27  %  enough  ArcAngles,  the  generated  path  is  unrealistic  and  ... 

results  in  a 

28  %  much  longer  flight  duration  than  is  realistic. 

29  % 

30  %  NOTE  2:  Within  the  function,  v_f  provides  an  approximation  . . . 

of  the 

31  %  maximum  velocity  which  the  subject  will  achieve  at  any  time  . . . 

along  its 

32  %  path.  Maximum  velocities  will  vary,  but  will  always  stay  within 

33  %  reasonable  bounds  (in  all  test  cases:  below  10  Km/s) .  A  very  . . . 

unique 
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34 


combination  of  inputs  would  be  required  to  produce  a  ... 
velocity  of  this 

35  %  magnitude  or  greater,  and  it  may  not  be  possible. 

36  % 

37 

38  %%  Define  constants 

39 

40  %  Check  to  see  if  a  time  series  is  included.  If  it  isn't,  create 

a  new 

41  %  scenario.  If  it  is,  then  generate  positions  based  on  times 

42  if  isstruct(sv) 

43  %  sv  is  already  populated 

44  [  sv  ]  =  generate_data  (  sv  )  ; 

45  else 

46  %  sv  needs  to  be  populated 

47  [  sv  ]  =  generate.sv  ( )  ; 

48  end 

49  end 

50 

51  function  sv  =  generate_sv  ( ) 

52  %  Generate  the  data  needed  to  get  sv  positions  at  any  time 

53 

54  %  Unit  vector  normal  to  a  plane  which  is  tangent  to  Earth 

55  N  =  randn  (3,1); 

56  N  =  N./norm(N); 

57 

58  %  Radius  of  the  Earth 

59  r  =  6378100; 

60 

61  %  The  number  below  is  the  minimum  arc  angle  for  which  an  sv  would 

be  used 

62  %  to  hit  a  subject.  The  full  definition  of  ArcAngle  is  given  as  a 

comment 

63  %  next  to  the  used  definition. 

64  ArcAngle  =  .  8 62325770 997 633+2 . 27 92 668825 921 60*rand;  %ArcAngle  =  . 

5500000/ r+(pi-5500000/r) *rand 

65 

66  Apogee  =  900000+rand*600000; 

67 

68  %%  Define  Ellipse 

69 

70  %  Y  coordiante  of  the  center  of  the  ellipse 

71  C_e  =  r-2800000; 

72 

73  %  The  constant  which  scales  the  y— components  of  the  ellipse 

74  a  =  Apogee  +  r  -  C_e; 

75 

76  %  This  determines  the  starting  point  of  the  subject's  path 

77  Beta  =  asin  (  (  C_e-r*cos  (ArcAngle/2  )  )  /a)  ; 

78 

79  %  The  constant  which  scales  the  x— components  of  the  ellipse 
so  b  =  r*sin  (ArcAngle/2  ) /cos  (Beta)  ; 
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81 

82  %%  Reentry  Phase 

83 

84  %  The  final  altitude  of  the  sv  during  the  reentry  phase 

85  a_f  =  150000+300000*rand; 

86 

87  %  The  velocity  at  the  end  of  the  reentry  phase 

88  v_f  =  5000+2000*rand; 

89 

90  Alpha_l  =  ... 

a  sin  (  (-a*  C_e  +  sqrt  (  (a*C_e)  ~  2-  (a" 2-b" 2 )  *  ( C_e  ~  2+b~2-  (a_f+r)  ~ 2 )  )  ) 

91  /  (a "2 -Jo' 2)  )  ; 

92  Alpha_2  =  ... 

a  sin  (  (-a*C_e-sqrt  (  (a*  C_e )  "2-  (a" 2-b" 2 )  *  ( C_e " 2+b"2-  (a_f+r)  '2)  )  ) 

93  /  (a"2-b"2)  )  ; 

94 

95  %  Ensure  that  the  correct  value  of  alpha  is  chosen  given  the  two 

possible 

96  %  values. 

97  if  (isreal (Alpha_l )  ==  0)  &&  ( isreal (Alpha_2 )  ==  0) 

98  error  (' Illegitimate  sv  Generation') 

99  elseif  isreal  (Alpha_l )  &&  isreal  (Alpha_2  ) 

100  if  alpha_l  <  alpha_2 

101  Alpha  =  Alpha_l; 

102  else 

103  Alpha  =  Alpha_2; 

104  end 

105  elseif  isreal  (Alpha_l )  ==  0 

106  Alpha  =  Alpha_2 ; 

107  else 

108  Alpha  =  Alpha_l ; 

109  end 
no 

in  omega  =  v_f  /  (sqrt ( (b*sin (Alpha) ) ~2+ (a*cos (Alpha) ) ~ 2)  *  2); 

112 

113  %  The  length  of  time  it  takes  to  complete  the  reentry  phase 

114  t_reentry_max  =  (Alpha+Beta) /omega; 

115 

116  %  Time  samples  across  the  orbital— phase 

in  t_max  =  (pi+2*Beta) / (2*omega) +omega*t_reentry_max; 

118 

119  sv.t_max  =  t_max; 

120  sv .  t_reentry_max  =  t_reentry_max  ; 

121  sv.a  =  a; 

122  sv.b  =  b; 

123  sv.  Alpha  =  Alpha; 

124  sv.Beta  =  Beta; 

125  sv.C_e  =  C_e; 

126  sv.  omega  =omega; 

127  sv.Rot_X  =  atan2  (N  (3)  ,  norm  (N  ( 1 :  2  )  )  )  ; 

128  sv.Rot_Y  =  rand*pi; 

129  sv.Rot_Z  =  — atan2  (N  (1)  ,N  (2)  )  ; 
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130  sv.N  =  N; 
m 

132  end 

133 

134  function  [  sv]  =  generate.data  (sv) 

135  %  Use  the  values  in  sv  to  generate  positions  &  velocities  at  .  . 

times  given  by  sv.t 

136 

137  reentry.length  =  sum(sv.t  <  sv.t_reentry.max); 

138 

139  t  =  sv . t ( 1 : reentry.length) ; 

mo  sv_r  =  — sv . r . *cos ( (t ./ sv . t_reentry_max) . *sv . omega . *t— sv . Beta) ; 
wi  sv_rv  =  2  *  sv.r  *  sv. omega  .*  t  .*  sin ( (t . /sv . t_reentry_max) . * 

142  sv  .  omega  .  *t— sv .  Beta)  . /sv . t.reentry.max; 

143  sv.r  (2,:)  =  sv.a.*sin(  (t./sv.  t .reentry  .max)  .  *  .  .  . 

144  sv .  omega  .  *t— sv .  Beta)  +sv  .  C_e ; 

145  s v_rv  ( 2 ,  : )  =  2  *  sv .  a  *  sv .  omega  .  *  t  .  *  .  .  . 

146  cos  (  (t ./ sv .  t_reentry_max)  .  *sv .  omega .  *t— sv .  Beta)  ... 

. /sv . t.reentry.max; 

147  %%  Orbital  Phase 

148 

149  t  =  sv .  t  ( reentry_length  +  l:end); 

150  sv_o  =  —  sv .  b  .  *cos  (2  .  *sv .  omega .  *t— sv .  Beta— sv .  omega*  ..  . 

151  sv  .  t_reentry_max)  ; 

152  sv.ov  =  2 *sv . b*sv .  omega  .  *sin  (2  .  *sv .  omega  .  *t— sv .  Beta—  ..  . 

153  sv  .  omega*  sv  .  t_reentry_max)  ; 

154  sv_o  (2,  : )  =  sv .  a  .  *sin  (2  .  *sv .  omega  .  *t— sv .  Beta-sv .  omega*  .  .  . 

155  sv  .  t_reentry_max) +sv  .  C.e  ; 

156  sv_ov(2,:)  =  2*sv .  a*sv .  omega  .  *cos  (2  .  *sv .  omega  .  *t-sv .  Beta  ..  . 

157  -sv .  omega*  sv  .  t_reentry_max)  ; 

158 

159  %%  Add  Phases 

160 

161  svp  =  sv.r  (1,  :)  ; 

162  svp  (2  ,  1 :  size  (  sv.r  ,  2  )  )  =  sv.r  (2,  :  )  ; 

163  svp  (1,  (size (sv.r,  2) +1)  :  (  size  (  sv.r  ,  2  )  +  size  (  sv.o  ,  2  )  )  )  ... 

164  =  sv.o  ( 1 ,  :  )  ; 

165  svp  (2,  (size(sv_r,2)+l)  :  (  size  (  sv.r  ,  2  )  +  size  (  sv.o  ,  2  )  )  )  ... 

166  =  sv.o  (2,  :  )  ; 

167  svp  (3,:)  =  zeros  ( 1 ,  size  ( svp,  2  ))  ; 

168 

i69  svv  =  sv.rv  ( 1 ,  :  )  ; 

no  svv  ( 2  ,  1 :  size  (  sv.rv ,  2  )  )  =  sv_rv(2,  :)  ; 

m  svv (1,  (size  (sv.rv, 2) +1)  :  (size  (sv.rv, 2) +size  (sv.ov, 2) ) )  .  .  . 

172  =  sv.ov  ( 1 ,  :  )  ; 

173  svv  (2,  (size  (sv.rv,  2)  +1)  :  (size  (sv.rv,  2)  +size  (sv.ov,  2)  )  )  .  .  . 

174  =  sv.ov  (2,  :  )  ; 

175  svv  (3,:)  =  zeros  ( 1 ,  size  (  svv,  2  ))  ; 

176 

177  sv.r  =  ZXY.Rotat  ion  ( svp,  sv  .  Rot.Z  ,  sv  .  Rot.X ,  sv  .  Rot.Y)  ; 
ns  sv.v  =  ZXY.Rotat ion ( svv,  sv . Rot.Z  ,  sv . Rot.X ,  sv . Rot.Y)  ; 

179  end 
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A.2  OrbGen.m 


1  function  [  Sat  ]  =  OrbGen (  N,  t  ) 

2  %%  Authorship 

3  %  Stephen  Hartzell,  AFIT  -  GRA 

4  %  stephen.hartzell.ctr@afit.edu,  hartzell.stephen@gmail.com 

5  %  Created:  02-04-2011,  Modified:  06—06—2012 

6  % 

7  %  OrbGen  Generate  satellite  orbital  paths. 

8  %  [  SAT  ]  =  OrbGen (  N,  T  )  generates  a  psuedo  random  realistic  . . . 

satellite 

9  %  orbit  on  the  opposite  side  of  an  Earth— tangental  plane  . . . 

defined  by  N. 

10  %  N  is  a  unit  vector  which  points  from  the  center  of  mass  of  ... 

the  Earth 

11  %  towards  a  plane  which  is  tangential  to  the  Earth.  T  is  an  array 

12  %  containing  times  at  which  satellites  positions  are  to  be  . . . 

generated . 

13  %  Positions  saved  to  SAT  correspond  to  the  times  in  T  are  . . . 

stored  as 

14  %  column  vectors  in  a  matrix  with  the  format  [Sx;  Sy;  Sz] . 

15  % 

16  %  Note:  T  is  restricted  less  it  become  impossible  to  generate  . . . 

an  orbital 

n  %  section  only  on  the  side  of  the  plane  opposite  the  Earth.  . . . 

For  this 

is  %  reason,  LEO  satellites  are  not  generated,  since  their  . . . 
orbital  periods 

19  %  are  comparable  to  the  flight  length  of  the  space  vehicle. 

20 

21  %%  Generate  Coherent  Sections  of  Orbits  Defined  By  Angular  Region 

22 

23  %  The  check  is  here  to  guard  against  the  rare  (and  perhaps 

24  %  impossible)  situation  of  orbit  parameters  along  with  N  causing 

25  %  the  function  to  be  unable  to  find  a  realistic  satellite  section 

26  chk  =  0  ; 

27  while  chk  ==  0 

28  switch  randi([2,4]) 

29  case  1 

30  %  LEO 

31  LEO.  name  =  'LEO'; 

32  LEO .  eccentricity  =  . l*rand; 

33  LEO .  SemiMa  jAx  =  7378000;  %  m 

34  LEO  .  SemiMinAx  =  LEO  .  SemiMa  jAx  *  ... 

sqrt ( 1-LEO . eccentricity ~ 2 ) ;  %  m 

35  LEO.inc  =  rand  *  90  *  pi/180; 

36  LEO .  LoAN  =  rand  *  360  *  pi/180; 

37  LEO.AoP  =  rand  *  90  *  pi/ 180; 

38  LEO.Vt  =  7800;  %  m/s 

39  LEO . Circumference  =  pi* (LEO . SemiMa jAx+LEO . SemiMinAx) .. . 

40  *(1  +  (  (3*  (  (LEO.  SemiMa  jAx-LEO.  SemiMinAx) /..  . 
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41 


42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 


(LEO . SemiMa jAx+LEO . SemiMinAx) ) "2 ) / ( 1 0  +  sqrt(4— 3*.. 
( (LEO. SemiMa jAx-LEO. SemiMinAx) /.  .  . 

(LEO. SemiMa jAx+LEO. SemiMinAx) ) "2) ) ) )  ; 

LEO.TrueAn  =  rand  *  360  *  pi/180; 

LEO. Period  =  LEO . Circumference  /  LEO.Vt;  %  s 
I  =  LEO; 
case  2 

%  MEO 

MEO . name  =  ' MEO ' ; 

MEO . eccentricity  =  .05*rand; 

MEO . SemiMa jAx  =  20200000;  %  m 

MEO . SemiMinAx  =  MEO . SemiMa jAx  *  ... 

sqrt ( 1— MEO . eccentricity ~ 2 ) ;  %  m 
MEO.inc  =  rand  *  90  *  pi/180; 

MEO . LoAN  =  rand  *  360  *  pi/180; 

MEO.AoP  =  rand  *  90  *  pi/180; 

MEO.Vt  =  4000;  %  m/s 
MEO . Circumference  =  ... 

pi* (MEO . SemiMa jAx+MEO . SemiMinAx) * ( 1  + . . . 

( ( 3* ( (MEO . SemiMa jAx— MEO . SemiMinAx) / (MEO . SemiMa jAx+ . 
MEO . SemiMinAx) )" 2 )/( 10  +  sqrt ( 4—3* ( (MEO . SemiMa jAx— . 
MEO . SemiMinAx) / (MEO . SemiMa jAx+MEO . SemiMinAx) ) "2 ) ) ) ) 
MEO.TrueAn  =  rand  *  360  *  pi/180; 

MEO. Period  =  MEO . Circumference  /  MEO.Vt;  %  s 
I  =  MEO; 
case  3 

%  GEO 

GEO . name  =  ' GEO ' ; 

GEO . eccentricity  =  . l*rand; 

GEO. SemiMa jAx  =  42164000;  %  m 

GEO . SemiMinAx  =  GEO . SemiMa jAx  *  ... 

sqrt ( 1-GEO . eccentricity " 2 ) ;  %  m 
GEO.inc  =  rand  *  90  *  pi/180; 

GEO. LoAN  =  rand  *  360  *  pi/180; 

GEO.AoP  =  rand  *  90  *  pi/180; 

GEO.Vt  =  3075;  %  m/s 
GEO . Circumference  =  ... 

pi* (GEO. SemiMa jAx+GEO . SemiMinAx) * (1  + . . . 

( (3* ( (GEO. SemiMa jAx-GEO. SemiMinAx) / (GEO . SemiMa jAx+ . 
GEO. SemiMinAx) ) "2) / (10  +  sqrt ( 4—3* ( (GEO . SemiMa jAx— . 
GEO. SemiMinAx) / (GEO . SemiMa jAx+GEO . SemiMinAx) ) ~2) ) ) ) 
GEO.TrueAn  =  rand  *  360  *  pi/180; 

GEO. Period  =  GEO . Circumference  /  GEO.Vt;  %  s 
I  =  GEO; 
case  4 

%  HEO 

HEO . name  =  ' HEO  ; 

HEO . eccentricity  =  .  2*rand; 

HEO. SemiMa jAx  =  60000000;  %  m 

HEO . SemiMinAx  =  HEO . SemiMa jAx  *  ... 

sqrt ( 1-HEO . eccentricity " 2 ) ;  %  m 
HEO.inc  =  rand  *  90  *  pi/180; 
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88 

HEO.LoAN  =  rand  *  360  *  pi/180; 

89 

HEO.AoP  =  rand  *  90  *  pi/180; 

90 

HEO.Vt  =  2600;  %  m/s 

91 

HEO . Circumference  =  ... 

pi* (HEO. SemiMa jAx+HEO . SemiMinAx) * (1  + . . . 

92 

( (3* ( (HEO. SemiMa jAx-HEO. SemiMinAx) / (HEO . SemiMa jAx+ . . . 

93 

HEO. SemiMinAx) ) ~ 2) / (10  +  sqrt (4—3* ( (HEO. SemiMa jAx— . . . 

94 

HEO . SemiMinAx) / (HEO . SemiMa jAx+HEO . SemiMinAx) ) ~ 

2))))} 

95 

HEO.TrueAn  =  rand  *  360  *  pi/180; 

96 

HEO. Period  =  HEO . Circumference  /  HEO.Vt;  %  s 

97 

I  =  HEO; 

98 

end 

99 

100 

I.ArcAng  =  (2*pi*max (t) ) /I .Period; 

101 

102 

%  Define  angle  between  V  and  N,  V  defines  the  orientation 

of  the 

103 

%  y— axis  of  the  plane  containing  the  generated  ellipse 

104 

if  strcmp ( I . name,  ' LEO ' ) 

105 

ang  =  rand*15*pi/180; 

106 

elseif  strcmp ( I . name, 'MEO') 

107 

ang  =  rand*45*pi/180; 

108 

elseif  strcmp ( I . name, 'GEO') 

109 

ang  =  rand*65*pi/180; 

110 

elseif  strcmp ( I . name, 'HEO') 

111 

ang  =  rand*65*pi/180; 

112 

end 

113 

114 

%  Generate  random  vector  U  and  find  the  the  cross  product 

of  ... 

N  and 

115 

%  U  to  create  a  unit  vector  which  over— writes  U  which  is  . 

.  . 

orthogonal 

116 

%  to  N. 

117 

U  =  randn (3,1) ; 

118 

U  =  cross (U, N) ; 

119 

U  =  U  . /  norm (U) ; 

120 

121 

%  Create  vector  V  which  defines  the  orientation  of  the  elliptical 

122 

%  orbits 

123 

U  =  6378000 . *tan (ang) . *  U; 

124 

V  =  6378000 . *N+U ; 

125 

126 

%  The  angles  which  define  the  orientation  of  the  elliptical  . . . 

orbits 

127 

I. LoAN  =  -atan2 (V ( 1 ) , V ( 2 ) ) ; 

128 

I . inc  =  atan2 (V ( 3 ) , norm (V ( 1 : 2 ) ) ) ; 

129 

I.AoP  =  rand*pi; 

130 

131 

%%  Find  Angular  Limits  for  Arguments  of  the  Ellipse's  Defintion 

132 

133 

%  Initialize  values 

134 

alpha  =  0; 

135 

A_alpha  =  .2; 
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136 

137 

138 

139 

140 

141 

142 

143 

144 

145 

146 

147 

148 

149 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 
161 
162 

163 

164 

165 

166 

167 

168 

169 

170 

171 

172 

173 

174 

175 

176 

177 

178 

179 

180 
181 
182 

183 

184 

185 

186 
187 


dist  =  2; 

%  This  finds  the  positive  angle,  alpha,  which  if  input  into  the 
%  equation  of  the  ellipse  in  its  own  plane  corresponds  to  one  of 
%  the  two  points  at  which  the  ellipse  intersects  the  equation  of 
%  the  plane  which  is  tangent  to  the  Earth  (as  defined  by  N) .  This 
%  loop  finds  alpha  by  an  iterative  approach, 
while  dist  >  le— 6 

%  Try  a  new  value  of  alpha 
alpha  =  alpha  +  A_alpha; 

%  Define  position  of  ellipse  in  2D  plane  at  alpha.  The 
%  convention  of  defining  the  ellipse  as  shown  below: 

%  x— values  =  —sin  and  y— values  =  —cos,  was  chosen  to  ensure 
%  that  for  alpha  =  0,  P_r  is  the  point  on  the  ellipse 
%  opposite  V  (see  V  above) 

Px  =  — I . SemiMa jAx*sin (alpha) ; 

Py  =  — I . SemiMinAx*cos (alpha) ; 

%  Rotate  the  previously  defined  positions  along  the  ellipse 
%  in  its  own  plane  to  define  the  points  P_r  in  3D  space 
%  which  correspond  to  the  points  in  the  ellipses  own  plane 
[  P_r  ]  =  ZXY_Rotation  (  [Px;Py;0],  I.LoAN,  I.inc,  I.AoP  ); 

%  ang  is  overwritten  and  now  contains  the  angle  between  the 
%  vector  pointing  from  the  center  of  mass  of  the  Earth  to  P_r 
%  and  N 

ang  =  acos ( (P_r . /norm ( P_r ) ) ' *N) ; 

if  ang  >  pi/2 

last_alpha  =  alpha; 
elseif  ang  <  pi/2 

A.alpha  =  A_alpha/2; 

%  Saves  the  last  angle  which  did  not  go  on  the  other 
%  side  of  the  plane 
alpha  =  last.alpha; 

end 

dist  =  abs (ang-pi/2 ) ; 

end 

%  Initialize  values 
beta  =  0; 
a  .beta  =  -.2; 
dist  =  2; 

%  This  finds  the  negative  angle,  beta,  which  if  input  into  the 
%  equation  of  the  ellipse  in  its  own  plane  corresponds  to  one 
%  of  the  two  points  at  which  the  ellipse  intersects  the 
%  equation  of  the  plane  which  is  tangent  to  the  Earth  (as 
%  defined  by  N) .  This  loop  finds  beta  by  an  iterative  approach. 
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188 

189 

190 

191 

192 

193 

194 

195 

196 

197 

198 

199 

200 
201 
202 

203 

204 

205 

206 

207 

208 

209 

210 
211 
212 

213 

214 

215 

216 

217 

218 

219 

220 
221 
222 

223 

224 

225 

226 

227 

228 

229 

230 

231 

232 

233 

234 

235 

236 

237 

238 

239 


while  dist  >  le— 6 

%  Try  a  new  value  of  beta 
beta  =  beta  +  A_beta; 

%  Define  position  of  ellipse  in  2D  plane  at  alpha.  The 
%  convention  of  defining  the  ellipse  as  shown  below: 

%  x— values  =  —sin  and  y— values  =  -cos,  was  chosen  to  ensure 
%  that  for  alpha  =  0,  P_r  is  the  point  on  the  ellipse 
%  opposite  V  (see  V  above) 

Px  =  — I . SemiMa jAx*sin (beta) ; 

Py  =  — I . SemiMinAx*cos (beta) ; 

%  Rotate  the  previously  defined  positions  along  the  ellipse 
%  in  its  own  plane  to  define  the  points  P_r  in  3D  space 
%  which  correspond  to  the  points  in  the  ellipses  own  plane 
[  P_r  ]  =  ZXY_Rotation  (  [Px;Py;0],  I.LoAN,  I.inc,  I.AoP  ); 

%  ang  is  overwritten  and  now  contains  the  angle  between  the 
%  vector  pointing  from  the  center  of  mass  of  the  Earth  to  P_r 
%  and  N 

ang  =  acos ( ( P_r . /norm ( P_r ) )  f  *N) ; 

if  ang  >  pi/2 

lastibeta  =  beta; 
elseif  ang  <  pi/2 

A.beta  =  A_beta/2; 

%  Saves  the  last  angle  which  did  not  go  on  the  other 
%  side  of  the  plane 
beta  =  last_beta; 

end 

dist  =  abs (ang-pi/2 ) ; 

end 

%  The  angular  region  in  which  the  satellite's  orbit  can  exist 
max.ang  =  2*pi  —  (alpha— beta) ; 

if  max.ang  <  I.ArcAng 
chk  =  0 ; 

else 

chk  =  1 ; 

%  Define  TrueAn  such  that  the  orbit  section  will  exist  on 
%  the  side  of  the  plane  pointed  to  by  N  and  is  randomly 
%  positioned  on  that  side 

I. TrueAn  =  alpha  +  rand  *  (max_ang  -  I.ArcAng); 

%  Angular  velocity 
Aw  =  2*pi  /  I. Period; 

%  Define  the  motion  of  the  satellite  in  its  own  plane 
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242 
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248 
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250 

251 
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Px  =  — I . SemiMa jAx*sin (Aw . *t  +  I.TrueAn); 

Py  =  — I . SemiMinAx*cos (Aw . *t  +  I.TrueAn); 

Pz  =  zeros ( 1 , length (Px ( 1 ,:))) ; 

%  Rotate  the  previously  defined  positions  along  the  ellipse 
%  in  its  own  plane  to  define  the  points  P_r  in  3D  space 
%  which  correspond  to  the  points  in  the  ellipses  own  plane 
[  P_r  ]  =  ZXY_Rotation  (  [Px;Py;Pz],  I.LoAN,  I.inc,  I.AoP  ); 

Sat(l,  (1 : size (P_r , 2) ) )  =  P_r(l,:);  %#ok<*AGROW> 
Sat(2,(l:size(P_r,2)))  =  P_r(2,:); 

Sat(3,  (l:size(P_r,2)))  =  P_r  (3,  : )  ; 

Sat  (4 , 1 : length (t ) )  =  t; 

end 

end 

end 


no 


Appendix  B:  Frames 


As  discussed  in  Section  2.4,  the  NLO  may  converge  to  any  one  of  several  solutions  to 
the  nonlinear  system  of  equations.  The  NLO  relies  on  the  assumption  that  the  measurement 
space  as  a  function  of  the  subject’s  state  £1  (x)  is  approximately  linear  near  the  initial  guess 
at  the  solution.  Early  in  this  research,  it  was  thought  that  this  linear  approximation  would 
be  more  accurate  if  the  angular  measurements  were  within  [-f,  f].  This  reasoning  is  no 
longer  considered  valid.  However,  for  this  reason,  the  NLOs  were  written  to  use  four 
different  coordinate  systems. 

It  was  thought  to  be  desirable  to  make  the  fl  (x)  function  approximately  linear.  From 
Equation  (3.9)  it  can  be  seen  that  the  LOS  measurements  depend  on  an  arctangent  function. 
From  Figure  B.l  it  is  evident  that  the  arctangent  function  is  well  approximated  by  a  line 
with  a  slope  of  one  on  the  interval  J-|,  |j.  However,  the  linear  approximation  may  be 
inaccurate  in  the  wings. 

Therefore,  a  set  of  four  coordinate  systems  are  used  to  keep  the  measurement  angles 
within  [-|,  |],  These  coordinate  systems  are  called  frames. 

If  the  LOS  measurements  fall  outside  of  the  desirable  interval,  then  the  measurements 
are  converted  to  another  frame,  and  the  optimization  operates  within  that  frame.  The  frames 
are  defined  as  shown  below. 

The  frames  in  which  the  measurements  are  placed  are  determined  as  follows.  The 
measurements  are  converted  to  frame 

1.  if  -f  <  0  <  f  and 

2.  if  -f  <  0  <  |  and  ^<6<^foxzjL<6<zf 

3.  if  |0|  >  |  and  -|<0<|or^<#<^ 
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Arctangent  Approximated  by  Line 


Figure  B.l:  Linear  approximation  of  arctangent  function.  The  arctangent  function  is 
reasonably  approximated  by  a  line  with  a  slope  of  one  over  the  domain  |]. 


4.  if  \(p\  >  |  and 

All  data  are  synthesized  in  frame  3.  When  using  the  NLO,  each  measurement  is 
converted  to  the  appropriate  frame. 
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Frame  1 
z-axis 


Frame  2 
z-axis 


y-axis 

x-axis 


y-axis 


Frame  3 
z-axis 


y-axis 


Frame  4 


z-axis 


y-axis 


Figure  B.2:  Frame  definitions.  The  four  frames  used  for  measurement  data  when 
performing  the  NLO. 
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Appendix  C:  Extended  Kalman  Filter 


This  appendix  represents  additional  research  that  was  conducted  to  improve  geoloca¬ 
tion  estimates.  This  approach  was  determined  to  be  inappropriate  because  it  assume  that 
estimates  of  the  subject’s  state  are  independent.  A  LOS  measurement  may  be  used  for  more 
than  one  estimation.  Therefore,  subject  estimates  are  not  independent.  However,  there  may 
still  be  some  utility  to  this  investigation,  and  so  it  is  included  here  in  the  case  of  future 
inquiry. 

Sensors  may  take  measurements  of  the  subject  in  question  over  long  periods  of  time. 
Therefore,  when  multiple  sensors  have  taken  measurements  at  nearly  the  same  time,  an 
estimate  of  the  subject’s  position  may  be  produced.  Therefore,  the  subject’s  track  will  be 
sampled  by  geolocation  estimates.  This  is  shown  in  Figure  C.  1 .  The  TA  and  NLO  considers 
each  geolocation  estimate  independently.  Therefore,  the  subject’s  track  is  never  taken  into 
consideration  by  the  NLO.  This  is  information  could  be  used  to  improve  the  geolocation 
estimates  and  confidence  estimates. 


Figure  C.l:  A  track  of  a  geolocation  estimates  each  with  their  own  covariance  matrices. 
There  is  an  apparent  track;  however,  neither  the  TA  nor  the  NLO  natively  take  this  track 
into  account. 


If  an  unknown  subject  is  being  observed,  it’s  track  is  likely  unknown;  however,  over 
short  lengths  of  time,  a  subject’s  track  may  be  accurately  described  by  a  second  order 
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model  where  acceleration  is  constant.  Uncertainty  in  this  model  may  be  described  by  its 
own  covariance  matrix.  Therefore,  an  geolocation  estimate  may  be  projected  forward  or 
backward  by  a  model  to  estimate  the  subject  at  another  time.  The  covariance  matrix  for  this 
new  projected  estimate  depends  on  the  confidence  in  the  original  estimate  and  the  model. 
The  Kalman  Filter  uses  this  scheme  to  improve  geolocation  estimates  at  all  points  in  time 
[22].  First  this  section  will  discuss  how  to  predict  and  refine  and  a  geolocation  estimate 
using  a  previous  estimate.  The  discussion  will  then  be  expanded  to  also  retrodict  from  the 
subsequent  estimate.  This  development  is  similar  to  that  given  in  [22]. 

C.l  The  Linear  Estimation  Problem 

The  linear  estimation  problem  involves  attempting  to  produce  the  ideal  estimate  of  a 
subject’s  state  given  an  imperfect  model  of  the  subject’s  motion  and  noisy  measurements 
of  this  subject’s  state.  The  subject’s  motion  is  described  as 

x/t+i  =  ^a-Xa  +  (C.l) 

In  this  equation,  <D/f  is  the  called  the  transition  matrix.  The  right  arrow  represents  that 
it  transforms  the  subject’s  state  forward  to  k  +  1.  It  describes  the  subject’s  motion  if  the 
model  were  perfect.  G*  describes  the  effect  of  noise  in  the  model  on  the  subject’s  next 
state.  For  our  purposes,  Gk  is  the  identity  matrix.  U*  is  a  vector  of  zero-mean  noise  with 
covariance  matrix  E®.  k  is  an  indexing  term. 

The  measurements  flk  associated  with  the  subject  at  xk  are  described  as 

n  =  M*x*  +  N,  (C.2) 

In  this  equation,  M*  is  a  matrix  that  transforms  the  subject’s  state  into  measurements. 
N*  is  a  noise  vector  with  covariance  matrix  In.  This  model  will  be  used  to  predict  the  state 
and  variance  at  k  +  1. 
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C.2  Prediction 


The  best  estimate  of  the  true  state  at  k  is  xk  =  xk+ek  where  ek  is  the  error  in  the  estimate. 
The  error  is  mean- zero  and  has  a  covariance  matrix  of  £Xjt.  Since  U/f  is  mean-zero,  the  best 
prediction  ~X£+1  is  given  by 


xk+i  =  Okxk. 


(C.3) 


To  find  the  covariance  matrix  of  the  predicted  estimate,  the  above  expression  is 
expanded  as 


x*+i  =  Okxk 


-  (Xk  +  6t) 


~  ®kxk  +  &kek 


Solving  Equation  (C.l)  for  <t>kxk,  and  substituting  this  into  the  above  equation  produces 


x/t+i  -  xk+\  ~  +  Okek 


(C.4) 


The  covariance  matrix  may  be  found  from  this  expression. 

The  expression  for  the  covariance  matrix  of  the  prediction  is  given  in  the  next  section 
in  the  set  of  five  Kalman  Filter  equations.  In  the  Backward-Forward  Extended  Kalman 
Filter,  the  only  equations  that  differ  are  the  prediction  equations  developed  in  this  section. 
While  the  prediction  equations  change  in  the  Backward-Forward  case,  these  results  are  still 
relevant. 

C.3  Kalman  Filter  Equations 

The  Kalman  Filter  is  described  performed  by  the  following  set  of  five  equations. 


116 


xk+1  =  a>A 


£jt+i  -  +  G^L(j)Gj 


Ki-l.  i  —  St+iMj 


-Ar+l^+l 


|^M^+|  L^+iM^+i  +  En 
X/t+l  =  Xk+I  -  Kk+l  (Mj.+  |  XJ-+!  -  fifc+l) 
2-ir+l  =  2fe+ 1  -  Kk+iyik+iHk+i 


(C.5) 

(C.6) 

(C.7) 

(C.8) 

(C.9) 


These  equations  are  given  in  [22].  These  equations  will  be  used  and  modified  for  the 
Backward-Forward  Extended  Kalman  Filter.  This  section  describes  how  to  use  prediction 
to  improve  the  estimate  of  the  subject’s  state  (the  geolocation  estimate)  and  to  improve  the 
confidence  in  this  estimate.  In  the  next  section,  this  idea  will  be  extended  to  incorporate 
the  subject’s  future  states  retrodicted  to  the  current  state  being  estimated. 

C.4  Adding  Retrodiction 

Rather  than  only  predicting  forward  to  the  subject’s  next  state,  the  Forward-Backward 
Extended  Kalman  Filter  also  retrodicts  to  the  subject’s  previous  state.  Therefore,  the 
prediction  and  retrodiction  must  be  used  together  to  provide  a  better  estimate  of  the 
subject’s  state  and  covariance  matrix. 

The  best  estimate  of  the  subject’s  state  at  k  -  1  and  at  k  +  1  are  described  by 


ii  - 

-  N  (xhTxjJ 

(C.10) 

Xfe+l  ' 

-  N (xfc+i,Hxt+1) 

(C.ll) 
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From  Equation  C.5,  it  can  be  can  seen  that  the  best  prediction  and  retrodiction  to  the 


subject’s  state  at  k  is  given  as 


x*  =  ®*_ix*_i  (C.12) 

X*  =  ®k+A+l  (C.13) 

Now  he  ideal  means  of  combining  the  prediction  and  retrodiction  must  be  found  to 
provide  the  best  estimate  of  x*  using  the  prior  and  posterior  subject  states.  This  combined 
estimate  is  called  x  k.  The  ideal  combined  estimate  is  the  subject  state  at  k  that  makes 
the  prediction  and  retrodiction  most  likely.  This  becomes  a  likelihood  estimation  problem 

^ _ y 

[18].  Therefore  x  k  may  be  found  as  given  by 

Vjt  =  arg  maxXfxVi,  x"*+i|x*)  (C.14) 

x* 

Expanding  this  equation  requires  knowledge  of  the  distribution  of  "x*_i  and  ~x&+i- 
These  distributions  are  Gaussian  with  covariances  which  may  be  calculated  as  shown  in 
Equation  (C.6).  Positive  constants  are  irrelevant  to  the  arg  max  problem,  so  they  are  ignored 
in 


x  k 


arg  max<? 


(- -x*)T£^i (Hk-i -x*)j  (- \ (W  -xa-)T 2A,ji(TA-+i -xt)) 


(C.15) 


X* 


The  arg  max  is  unaffected  by  taking  the  natural  log  of  the  argument.  Therefore,  by 
taking  the  natural  log  of  this  argument  and  eliminating  the  factor  of  -1/2,  the  equation 
is  simplified  to 


x  k 


arg  min  p*_i  -  xt)  p*-i  -  x*)  +  (x~*+i  -  x*)  (x~t+i  -  x*)  (C.16) 

Xk 
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This  equation  may  now  by  solved  for  x  k  by  finding  the  value  of  xk  such  that  the  first 
derivative  of  the  argument  is  zero.  Calculating  the  derivative  requires 

Vx(xTMx)  =  (MT  +  M)x  (C.17) 

Therefore,  by  using  Equation  C.17  on  Equation  C.16  and  setting  the  result  to  zero, 
yielding 


+£?-,)(*-,  "**)  +  (£w  +  - x»)  =  0  (C.18) 

Because  covariance  matrices  are  symmetric  by  their  definition,  this  result  may  be 
further  simplified  to 


2 (xx-i  -  x*)  +  (Vv+1  -  xfc)  =  0 

£*- “  S*- ixa  +  ~  ^a+ixa  =  0 

^A-1XA  +  E^X*  +  S^iX^i 

+  ^A+l  )XA  =  ^k- 1  X*-l  +  ^A+l  X-v+l 
XA  =  +  ^A+l)  Xx-1  +  ^A+l  x-r+l  j 

~a  =  (e^  +  E^)'  1  (e*-!  x*-!  +  E^twl)  (C.19) 

The  last  step  follows  from  Equation  C.14.  This  equation  provides  a  means  of 

calculating  the  optimal  prediction  of  x*.  Using  this  result,  the  covariance  matrix  of  the 
^ ^ 

combined  prediction  E  k  can  be  calculated.  The  definition  of  the  covariance  matrix  is  used 
^ ^ 

to  find  E  *  as  given  by 


2  a 


(C.20) 
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The  k  -  E  [*x%]  term  will  be  tackled  first.  It  will  be  expanded  and  errors  in  the 

state  estimate  at  k  -  1  and  k  +  1  will  be  taken  in  consideration  as  in  Section  C.2.  Before 
^ ^ 

examining  this  term,  x  k  is  expanded  out  as  follows  via  Equation  C.5  and  Equation  C.19 
as  given  by 


Let  A  =  (2^1+ £^,) 

V,  =  A  (C.21) 

From  Equations  (C.3)  and  (C.4),  can  be  expanded  O^x*^  and  ^+1xt+1  as 
x  k  =  A (xk  -  G^._ i U,_ i  +  <5>k-\ek-\  j  +  ^£+1  (xfe  _  Gyt+iUfe+i  +  dE+i e,-+|  jj 

=  ^  ^jfc-lXi:  +  ^jt+lXttj  +  A  ~  Gjt-lUjt-l  j  +  -  Gfc+lUyt+l  jj 

(C.22) 

The  left-hand  term  in  the  final  equation  is  constant  since  xk  is  the  true  value  of  x  at 
k.  The  right-hand  term  mean  0  because  ek  |  and  ek  +  1  are  mean-zero.  Recall  that  it  has 
been  shown  in  Equation  (2.53)  that  the  NLO  is  an  unbiased  estimator.  Therefore,  in  the 
*x%  -  E  [^x^J  term,  the  constant  factors  are  subtracted  out  since  the  expected  value  of  a 
constants  yields  the  same  constant.  The  expected  value  of  the  right  hand  term  is  zero. 
Therefore, 

x^-E^x^j=A  |  -  Gjt-iUfc_i  j  +  'E‘k+\  ^/t+i^+i  _  Gfe+iUfe+i  jj .  (C.23) 

Substituting  this  result  into  Equation  (C.20)  produces  an  equation  too  large  to  write 
out  here.  Therefore,  the  result  of  this  substitution  will  be  shown  by  placing  the  U*_i  and 
Ujfc+i  terms  into  a  variable  H.  The  expansion  of  the  substitution  shown  below  is  used  to 
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demonstrate  that  each  term  in  Equation  (C.23)  may  be  analyzed  separately,  ie 


=  E 


=  E 


Xk  =  A  E 


A  +  Efc+i^+i^+i  +  rU]j\^y-‘k-\®k-iek-\  +  ^fc+i^t+iefc+i  +  ^ 

+«Wk_, +Ut)a 


’Ekll®k-\£k-i  +  S^jO^+ie^+i  +'W](6AT_,<Dj_1i:jti1  +  e^iO^X;^!  +  'W1 


(C.24) 


Due  to  the  linearity  of  expected  value  [7],  the  expected  values  of  the  cross-terms  may 
be  examined  separately.  Since  the  only  random  variables  present  are  e^_,  and  ek+l,  the 
expected  value  will  only  operate  on  these.  One  of  the  cross-terms  is 


E^i%_iE[e,_i6tT+1]Oj+1i: 


1-1 

*k+l 


(C.25) 


Note  that  from  [7].  The  covariance  matrix  of  two  random  vectors  X  and  Y  may  be 
described  as 


X(X,  Y)  =  E  [XYT]  -  E[X\E  [T]t  (C.26) 

If  X  and  Y  are  mean-zero,  then  this  simplifies  to 

E(X,  Y)  =  E  [XYT]  (C.27) 

since  the  subject’s  state  estimates  are  mean-zero  and  are  based  on  independent  identically 
distributed  measurements,  their  covariance  is  zero.  This  same  argument  holds  for  the 
expected  value  of  the  U  and  e  cross-terms.  Therefore,  Equation  (C.24)  simplifies  to  the 
expected  value  of  the  like-terms  as  shown  below. 
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^ ^ 

=  AE 


s  *-  i <t>k- 1 <T- 1 1 1 ^ ^ i  +  E^+1^+1Oj+1E^  ^CUCW 

=  A  [ei_.^_1]Oj_1E?i1  +  %>**,£  [e^*,]  S^i,  + 

=  A  ®:_,  „  0^,2^,  +  WWT)  AT  (C.28) 


Likewise,  expanding  out  and  simplifying  the  <U<UT  term  produces 


+  E^G^G^E^ 


(C.29) 


We  now  have  a  means  of  calculating  the  covariance  matrix  of  the  combined  prediction  and 
retrodiction  estimate  of  the  subject’s  state.  Equations  C.28  and  C.19  may  be  used  in  the 
Kalman  Filter  equations  given  in  Section  C.3. 

In  this  way,  the  subject’s  track  is  taken  advantage  of  to  refine  the  state  estimates  and 
the  confidences  in  those  estimates.  The  confidence  in  the  transition  matrix  may  not  be 
known.  Therefore,  it  is  common  to  use  a  guess  at  the  first  iteration  of  the  Kalman  Filter. 

This  approach  was  deemed  inappropriate  since  the  same  FOS  measurements  may  be 
used  to  geolocate  the  subject  at  different  times.  Therefore,  the  geolocation  estimates  are 
not  linearly  independent  thus  violating  a  Kalman  Filter  assumption. 
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